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PREFACE 


The  study  reported  herein  was  sponsored  by  the  Office,  Chief 
of  Engineers,  U.  S.  Army,  as  part  of  the  Civil  Works  General 
Investigations,  Environmental  and  Water  Quality  Operational 
Studies  ( EWQOS )  Program.  Work  Unit  No.  31593  (Task  IA.4)  enti¬ 
tled  "Improve  and  Verify  Multidimensional  Hydrodynamic  Mathemat¬ 
ical  Models  for' Reservoirs"  supported  the  subject  study. 

The  study  was  conducted  during  the  period  July  1978  to 
October  1980  by  Dr.  John  E.  Edinger  and  Mr.  Edward  M.  Buchak  of 
J.  E.  Edinger  Associates,  Inc.  Mr.  Mark  S.  Dortch  and  Dr.  Billy 
H.  Johnson  of  the  Hydraulics  Laboratory,  U.  S.  Army  Engineer 
Waterways  Experiment  Station  (WES),  monitored  the  effort.  This 
report  was  written  by  Mr.  Buchak  and  Dr.  Edinger.  Program  man¬ 
ager  of  EWQOS  was  Dr.  Jerome  L.  Mahloch,  WES  Environmental 
Laboratory . 

Commanders  and  Directors  of  WES  during  this  study  and  the 
preparation  of  this  report  were  COL  John  L.  Cannon,  CE,  and  COL 
Tilford  C.  Creel,  CE.  Technical  Director  was  Mr.  Fred  R.  Brown. 


This  report  should  be  cited  as  follows: 


Buchak,  E.  M. ,  and  J.  E.  Edinger.  1982.  "User  Guide  for 
LARM2 :  A  Longitudinal-Vertical,  Time-Varying  Hydrodynamic 
Reservoir  Model,"  Instruction  Report  E-82-3  ,  prepared 
by  J.  E.  Edinger  Associates,  Inc.,  for  the  U.  S.  Army 
Engineer  Waterways  Experiment  Station. 


Accession  For 

NTIS  GRAAI 
DTIC  TAB 


Unannounced  □ 

Justification _ 


By - 

Distribution/ 


Availability  Codes 


Dist 


Avail  and/or 
Special 


TABLE  OF  CONTENTS 


PREFACE  .  i 

LIST  OF  FIGURES . iii 

1.  INTRODUCTION .  1 

2.  CAPABILITIES  AND  LIMITATIONS .  3 

3.  THEORETICAL  BASIS  .  5 

4.  APPLICATION  PROCEDURE  AND  EXAMPLE  .  11 

4.1  Required  Data .  11 

4.2  Geometric  Schematization  .  12 

4.3  Time-Varying  Boundary  Condition  Data  .  17 

4.4  Initial  Conditions,  Hydraulic  Parameters,  and 

Constituent  Reactions .  20 

4.5  Input,  Output  and  Computer  Resource 

Requirements  .  25 

4.6  Dillon  Reservoir  Example  .  27 

REFERENCES .  32 

APPENDIX  A  LARM2  Input  Data  Description 
APPENDIX  B  Example  Input  Data 
APPENDIX  C  Example  Output  Data 
APPENDIX  D  LARM2  Error  Messages 


APPENDIX  E  LARM2  Glossary  of  Important  FORTRAN  Variables 

APPENDIX  F  Preprocessor  GIN 

APPENDIX  G  Postprocessor  VV 

APPENDIX  H  List  of  LARM  Application  Reports 

APPENDIX  I  Julian  Data  Calendar 


i 


LIST  OF  FIGURES 


FIGURE  1  Locations  and  Sign  Conventions  of 
Major  FORTRAN  Variables  on  Finite 
Difference  Grid . 

FIGURE  2  Dillon  Reservoir  Finite  Difference 

Grid . 

FIGURE  3  Dillon  Reservoir  . 

FIGURE  4  Dillon  Reservoir  LARM  Segmentation  .  . 


1 


USER  GUIDE  FOR  LARM2 :  A  LONG I TUP INAL- VERT I CAL. 

TIME-VARYING  HYDRODYNAMIC  RESERVOIR  MODEL 

1.  INTRODUCTION 

This  report  is  the  user  guide  for  a  code  that  permits  time- 
varying  hydrodynamic  and  transport  simulations  of  rivers,  lakes, 
and  impoundments.  The  code  is  called  LARM2,  an  acronym  for 
laterally-averaged  reservoir  model,  version  two.  The  original 
LARM  was  developed  for  the  U.  S.  Army  Engineer  Division,  Ohio 
River.  Three  reports  describing  its  development  have  been  writ¬ 
ten,  the  final  one  of  which  summarizes  the  mathematical  basis  of 
the  code.  That  report,  A  Hydrodynamic,  Two-Dimensional  Reservoir 
Model:  Development  and  Test  Application  to  Sutton  Reservoir  Elk 
River,  West  Virginia  (Edinger  and  Buchak,  1979),  also  includes  a 
discussion  of  the  code's  first  application,  with  a  comparison  of 
model  results  and  prototype  data. 

Since  that  time,  a  dozen  or  more  applications  of  LARM  have 
been  made  and  additional  development  work  has  been  funded  by 
U.  S.  Army  Engineer  Waterways  Experiment  Station.  The  newer 
version  of  the  code,  LARM2,  incorporates  features  that  simplify 
its  use  and  expand  its  capabilities.  Among  these  are 

*  additional  coding  to  permit  transport  computations  for 
water  quality  constituents  beyond  temperature; 

*  the  ability  to  expand  and  contract  the  finite  difference 
grid  in  both  the  vertical  and  longitudinal  directions  as 
reservoir  volume  changes  in  response  to  flood  events  and 
drawdown ; 

*  a  decrease  of  30%  in  central  processor  time  requirements 


as  a  result  of  programming  improvements  to  accommodate 
vector-processing  central  processors; 

‘  reduction  of  application-related  code  changes  to  modify  a 
subroutine  for  the  input  of  time-varying  boundary  condi¬ 
tion  data  and  specify  constituent  reaction  rates; 

'  the  development  of  a  preprocessor  code,  GIN,  to  prepare 
LARM2  geometric  data  from  the  output  of  the  Hydrologic 
Engineering  Center  program  GEDA  (Hydrologic  Engineering 
Center,  1976);  and 

the  development  of  a  postprocessor  code,  VV,  to  plot 
velocity  vectors  using  LARM2  output. 

This  user  guide  contains  an  overall  view  of  LARM2  (Chapter 
2),  a  summary  of  its  theoretical  basis  (Chapter  3),  and  an  example 
problem  (Chapter  4).  A  thorough  understanding  of  the  example  as 
well  as  the  input  data  description  (Appendix  A)  is  required  of 
the  user  who  is  preparing  an  application. 


3 


2.  CAPABILITIES  AND  LIMITATIONS 

LARM2  was  developed  to  assist  in  the  analysis  of  water 
quality  problems  in  rivers,  lakes,  and  reservoirs  where  buoy¬ 
ancy  is  important  and  where  lateral  homogeneity  can  be  assumed. 

The  code  is  recommended  for  those  cases  where  longitudinal  and 
vertical  temperature  or  constituent  gradients  occur.  LARM2 
generates  time-varying  velocity,  temperature,  and  water  quality 
constituent  fields  and  surface  elevations  on  a  longitudinal  and 
vertical  grid.  Both  the  spatial  and  temporal  resolution  can  be 
varied  by  the  user,  within  certain  limits.  Units  used  are  the 
International  System  of  Units  (SI);  i.e. ,  metre-kilogram-second. 

An  application  requires  a  considerable  effort  in  planning 
and  data  assembly  and  computer  resources.  The  user  is  re¬ 
quired  to  modify  a  subroutine  that  supplies  time-varying  boundary 
condition  data  to  LARM2  for  his  particular  case.  If  the  user 
intends  to  make  use  of  the  constituent  transport  computation 
option,  he  must  be  able  to  quantify  the  reaction  rates  for  each 
constituent  in  terms  of  every  other  constituent  and  code  these 
statements.  The  user  must  also  be  able  to  review  results  for 
reasonableness  and  applicability  to  his  problem. 

LARM2  in  its  present  form  can  be  applied  to  a  single,  con¬ 
tinuous  reach  of  a  river,  lake,  or  reservoir.  Other  configura¬ 
tions  are  possible,  such  as  carrying  the  computations  into 
branches  and  tributaries,  if  the  code  is  employed  as  a  subroutine 


4 


for  each  reach  and  boundary  conditions  between  reaches  are 
specified.  This  configuration  has  not  been  tested,  however, 
and  would  require  significant  additional  programming. 

The  code  automatically  increases  the  number  of  horizontal 
layers  when  water  surface  elevations  increase  and  then  adds 
segments  as  backwaters  progress  into  the  upstream  end  of  a 
reservoir.  Similarly,  the  grid  contracts  both  vertically  and 
longitudinally  during  drawdown. 

Waterbody  simulations  through  long  periods  are  economically 
feasible.  An  eight -month  simulation  at  a  time  step  of  15  minutes 
for  temperature  and  one  other  constituent  on  a  grid  that  is 
30  segments  by  20  layers  with  375  active  cells  would  cost  $300 
on  a  commerical  batch-processing  computer  system.  Each  additional 
constituent  would  add  approximately  $30  to  this  cost.  These 
cost  estimates  are  exclusive  of  any  intermediate  simulations  and 
postprocessing  of  results.  Experience  has  shown  that  the  former 
may  account  for  ten  simulations  for  each  final  simulation  obtained. 

LARM2  has  a  simple  algorithm  for  the  computation  of  ice  for¬ 
mation  and  breakup  that  does  not  explicitly  account  for  the 
latent  heat  content  of  the  ice.  Typical  LARU2  simulations  will 
begin  in  March  and  continue  through  October  for  a  temperate  zone 
waterbody.  With  care  in  evaluating  computed  ice  cover,  these 
simulations  may  be  extended  to  the  November-to-February  period. 

LARM  and  LARM2  have  been  tested  against  analytic  solutions 
and  have  been  applied  to  a  dozen  or  more  field  cases.  A  list  of 
reports  describing  these  applications  is  given  in  Appendix  H. 


3.  THEORETICAL  BASIS 


Details  of  the  formulation  of  the  governing  partial  differ¬ 
ential  equations  and  the  subsequent  casting  of  those  equations 
into  finite  difference  form  are  presented  in  Edinger  and  Buchak, 
1979.  Briefly,  six  equations  (longitudinal  momentum,  vertical 
momentum,  continuity,  heat  and  constituent  balances,  and  state) 
are  solved  to  obtain  longitudinal  and  vertical  velocity  components, 
surface  elevations,  temperatures,  constituent  concentrations,  and 
densities  at  a  grid  of  points  in  space  and  time.  The  three- 
dimensional,  time- varying  partial  differential  equations  are  for¬ 
mally  averaged  over  the  reservoir  width  and  then  cast  into  finite 
difference  form  after  vertical  averaging  over  a  horizontal  layer 
thickness,  h,  with  boundaries  at  z=k+i  and  z=k-£.  The  laterally 
averaged  equations  as  vertically  integrated  over  the  layer  thick¬ 
ness  are  presented  below: 

longitudinal  (x-direction)  momentum 

w  <™h>  +  if  <D’Bh>  ♦  Ww  -  <"bV>k-i, 

*  k  s <™»  -  *,  £  wbw  +  -  °  <i*> 

with 


C*  p  /p  W  1  cos  $ 

&  & 

(surface) 

(lb) 

-A  3U/3z 
z 

(interlayer) 

(lc) 

g  u|u|/c2 

(bottom) 

(Id) 
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continuity 


<V»H  *  <V>k*  +  Si  ,UB«  -  H 


3 (Zb)  _  _9 
3t  3x 


(UB)dz  +  ||  -  0 


(internal) 
(over  total  depth) 


(2a) 

(2b) 


H 


vertical  (z-direction)  momentum 


3P 


31"  pg 


heat  balance 


si  <Bh«  +  s!  <«*“>  +  <»bbI)rt  -  <v»kJj 

_  3  f  3BhTx  _  (  3BT.  ,  3BT.  _  Hn  Bh 

3x  (Dx  3x  }  (Dz  3z  +  (Dz  ~ F 


constituent  balance 


si  <Bhc>  +  si  (DBhc>  +  <"bbc>l^  -  <»bbc>k-i, 


_i  rn  9BhC 
'  3x  ^  x  3x 


3BC. 


Hn  Bh 


.)  _  /D  zaz.-  +  (D  1®C\  ,  fin. 

;  lBz  3z  k-Us  V  z  3z  >k->s  V 


state 


p  -  5.29157  x  10-5  T3  -  8.45123  x  l(f3  TZ 


(3) 


(4) 


(5) 


I  -  -4 


■■  ■:  A 


+  6.59583  x  io“2  T  +  999.841 


(6) 
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where 

Ax  x-direction  momentum  dispersion  coefiicient  (m2/s) 

A_  z-direction  momentum  dispersion  coefficient  (m2/s) 

b  lake  width  (m) 

B  laterally  averaged  lake  width  integrated  over  h  (m) 

c  Chezy  resistance  coefficient,  m^/s 

C  laterally  averaged  constituent  concentration  inte¬ 

grated  over  h  (mg-1”1) 

C*  resistance  coefficient 

D  x-direction  temperature  and  constituent  dispersion 

x  coefficient  (m2/s) 

D  z-direction  temperature  and  constituent  dispersion 

z  coefficient  (m2/s) 

g  acceleration  due  to  gravity  (m/s2) 

h  horizontal  layer  thickness  (m) 

H  total  depth  (m) 

H  source  strength  for  heat  balance  (°C*m3  •  s-1 )  or  con¬ 

stituent  balance  (mg *1-1 • m3 • s-1 ) 

k  integer  layer  number,  positive  downward 

P  pressure  (Pa  =  N/m2 ) 

Q  tributary  inflow  and  withdrawal  rates  (m3/s) 

t  time  (s) 

T  laterally  averaged  temperature  integrated  over  h  (°C) 

U  x-direction,  laterally  averaged  velocity  integrated 

over  h  (m/s) 

ub  x-direction,  laterally  averaged  velocity  (m/s) 

V  cell  volume  (B'h-Ax)  (m3). 

Wa  wind  speed  (m/s) 

wb  z-direction,  laterally  averaged  velocity  (m/s) 
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x  and  z  Cartesian  coordinates:  x  is  along  the  lake  centerline 
at  the  water  surface,  positive  to  the  right,  and  z  is 
positive  downward  from  the  x-axis  (m) 

Z  surface  elevation  (m) 

p  density  (kg/m3) 

Pa  air  density  (kg/m3) 

t  z-direction  shear  stress  (m2/s2) 

z 

<j>  wind  direction  (rad) 

The  finite  difference  operation  is  applied  to  Equations 
1  through  5  and  introduces  two  additional  variables,  Ax,  the 
x-direction  spatial  step  (m),  and  At,  the  computation  time 
step  (s).  This  operation  produces  finite  difference  equation 
analogues  of  Equations  1  through  5.  The  solution  technique  is 
to  substitute  the  forward  time,  longitudinal  momentum  term  from 
Equation  la  into  the  vertically  integrated  continuity  expression, 
Equation  2b,  to  give  the  free  surface  frictionally  and  inertially 
damped  longwave  equation.  The  latter  is  solved  implicitly  for 
the  surface  elevation,  Z,  which  eliminates  the  longwave  speed 
time  step  restriction  (Ax/At  >  /g5) .  However,  the  computation 
time  step>At,  is  still  limited  by  the  fundamental  condition  that 

At  <  V/Q*  (7) 

for  each  cell  where  Q*  is  the  flow  into  or  out  of  the  cell  and 
V  is  the  cell  volume. 

The  computation  proceeds  as  follows.  New  water  surface 
elevations  are  computed  using  the  implicit  equation  combination, 
Equations  la  and  2b.  New  longitudinal  velocities  on  the  grid 


are  then  computed  from  the  explicit  Equation  la.  Equation  2a 
is  used  to  compute  vertical  velocities,  and  Equation  2b  is  used 
to  check  water  balances  prior  to  the  implicit  solution  of 
Equations  4  and  5  for  temperature  and  constituent  concentrations. 
Finally,  densities  are  computed  from  Equation  6*.  Results  may 
then  be  printed  and  a  new  time  step  computation  begins  once  again 
with  the  solution  for  water  surface  elevation. 

The  FORTRAN  names  of  the  variables  introduced  in  Equations 
1  through  5  are  given  below**: 


FORTRAN  name 
QTRIB,  QWD 
Tl,  T2 
U 

VOL 

WA 

W 

Zl,  Z2 

DLT 

DLX 

RHO 

RHOA 

ST 

PHI 


*  Personal  Communication,  August  1974,  Alan  Toblin,  NUS  Corpora 
tion,  Rockville,  Md. 

**  A  glossary  of  the  FORTRAN  variables  discussed  in  this  report 
is  presented  in  Appendix  E. 


The  location  of  the  variables  on  the  finite  difference  grid 
is  important  in  understanding  the  FORTRAN  coding  of  LARM2  and  the 
output.  Figure  1  shows  the  location  of  velocities  and  dispersion 
coefficients  at  cell  boundaries.  This  convention  permits  bound¬ 
ary  values  of  these  variables  to  be  exactly  zero.  Other  variables 
are  located  at  the  cell  center  and  represent  averages  for  the 
entire  cell. 

Throughout  the  computation  new  values  of  variables  replace 
old  values.  The  solution  technique  requires  that  values  of  Z, 

H,  T,  and  C  be  retained  for  two  time  steps.  At  the  end  of  a 
time  step,  the  contents  of  these  variable  arrays  are  exchanged 
so  that  the  older  values  are  placed  in  arrays  with  a  "1"  suffix 
and  the  newer  values  are  placed  in  arrays  with  an  "2"  suffix. 
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4.  APPLICATION  PROCEDURE  AND  EXAMPLE 

Applying  LARM2  to  a  particular  case  requires  five  steps: 

(1)  assembling  the  necessary  data;  (2)  schematizing  the  res¬ 
ervoir  geometry;  (3)  modifying  a  subroutine  for  the  input  of 
time-varying  boundary  condition  data;  (4)  setting  hydraulic 
and  other  parameters  and  constituent  reaction  rates;  and 
(5)  executing  the  code  and  evaluating  the  results.  These  steps 
are  described  below. 

4.1  Required  Data 

The  initial  task  in  applying  LARM2  to  a  specific  case 
(after  deciding  that  the  problem  requires  a  two-dimensional, 
time-varying  simulation  for  its  solution)  is  assembling  the 
necessary  data.  Reservoir  geometry  is  of  immediate  usefulness 
and  should  include  bathymetric  cross-sections  (y-z  plane). 

A  plan  (x-y  pl...:»e),  an  elevation  (x-z  plane),  and  an  elevation - 
area-volume  table  are  useful,  but  not  necessary,  adjunct 
information.  These  data  will  be  used  to  construct  the  compu¬ 
tational  grid. 

Initial  and  boundary  condition  data  are  required.  The 
former  consists  of  an  initial  water  surface  elevation  and  tem¬ 
perature  for  the  starting  day  of  the  simulation.  Boundary 
condition  data  in  their  most  general  form  consist  of  time-varying 
inflow  rates  and  temperatures  for  the  upstream  inflow  and  each 
significant  tributary,  time-varying  release  and  withdrawal  rates, 


and  time-varying  meteorological  data  for  wind  shear,  evaporation, 
and  surface  heat  exchange  calculations.  As  LARM2  is  now  coded, 
surface  heat  exchange  is  computed  using  equilibrium  temperature, 
the  coefficient  of  surface  heat  exchange,  and  solar  radiation. 
These  parameters  are  computed  separately  from  air  temperature; 
dew  point  temperature;  wind  speed;  and  observed  solar  radiation, 
percent  of  possible  sunshine,  or  cloud  cover  data  (Edinger,  Brady, 
and  Geyer  1974).  In  their  simplest  form,  the  boundary-condition 
data  are  constants.  For  each  water  quality  constituent  to  be 
computed,  initial  and  boundary  condition  data,  as  well  as  reaction 
rates,  are  needed.  Hydraulic  properties — Chezy  and  dispersion 
coefficients — and  solar  radiation  absorption  and  attenuation 
characteristics  are  also  required. 


4.2  Geometric  Schematization 

The  geometric  data  are  organized  into  a  grid  like  the  one 
shown  in  Figure  2  for  Dillon  Reservoir.  The  basic  parameters  to 
be  selected  in  defining  the  grid  are  the  longitudinal  spacing, 

Ax  (m),  and  the  vertical  spacing,  h  (m).  These  parameters  are 
constant  throughout  the  grid  and  define  two  of  the  three  di¬ 
mensions  of  each  cell.  The  third  dimension,  the  cell  widths, 
are  most  easily  taken  from  cross-section  drawings  and  represent 
an  average  value  for  each  cell.  The  vertical  columns  defined 
by  Ax  are  called  segments,  and  the  horizontal  rows  defined  by  h 
are  called  layers. 
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A  more  sophisticated  procedure  for  obtaining  cell  widths 
involves  the  use  of  the  Hydrologic  Engineering  Center  code  GEDA 
and  the  preprocessor  GIN,  the  user  notes  for  which  are  included 
as  Appendix  F  in  this  guide,  GEDA  utilizes  cross-sectional  data 
from  transects  at  irregular  intervals  and  produces  reservoir 
widths  as  a  function  of  elevation  at  regular  intervals.  GEDA 
also  computes  areas  and  volumes  as  functions  of  elevation. 
Cross-sectional  data  may  be  modified  using  GEDA  procedures  in 
order  to  fit  the  elevation-area-volume  table  to  one  computed 
from  detailed  topography  and  published  in  the  design  report  for 
the  reservoir.  GEDA  allows  several  Ax's  and  h's  to  be  tried 
with  only  small  additional  effort.  GIN  uses  GEDA  output  as  input 
and  produces  the  geometry  (BA)  cards  read  by  LARM2 .  GIN  also 
may  be  manipulated  to  change  the  computed  elevation-area-volume 
table  to  more  closely  match  the  given  table.  GIN  produces  an 
active/ inactive  cell  map  to  allow  screening  of  illegal  geometric 
features . 

Implicit  in  the  geometric  schematization  is  the  identification 
of  the  waterbody  center  line.  This  may  be  defined  as  the  dominant 
flow  path,  with  its  origin  at  the  upstream  boundary.  The  center 
line  may  follow  a  path  that  coincides  with  the  thalweg.  All  lon¬ 
gitudinal  distances  are  measured  along  the  selected  center  line. 

The  choice  of  Ax  and  h  is  important  to  the  success  of  the 
simulation.  By  selecting  the  ratio  h/Ax  to  match  the  overall 
reservoir  slope,  the  user  will  find  that  the  reservoir  bottom 
may  be  modelled  more  accurately.  While  relatively  small  values 


of  Ax  and  h  permit  finer  resolution,  the  resulting  large  number 
of  grid  cells  requires  more  computation  and  storage.  A  second 
consideration  in  the  choice  of  Ax  and  h  is  the  basic  stability 
criterion  that  At,  the  integration  time  step,  be  less  than  the 
volume  of  each  cell  divided  by  the  flow  through  the  cell 
(Equation  7).  This  criterion  also  affects  economy  by  forcing 
smaller  At's  as  cell  dimensions  are  decreased. 

In  practice,  At  can  be  changed  for  a  particular  simulation, 
while  the  geometry  established  initially  remains  fixed.  Sur¬ 
prises  from  an  implicit  choice  of  At  can  be  avoided  by  antic¬ 
ipating  the  limiting  At  dictated  by  the  geometry.  In  most 
cases,  this  limiting  At  will  occur  at  tributary  or  upstream 
inflows  and  withdrawals  or  downstream  outflows.  Equation  7  can 
be  evaluated  using  the  smallest  cell  volume  and  largest  antic¬ 
ipated  flow  at  these  inflow  or  outflow  segments. 

The  code  automatically  expands  the  grid  to  accommodate 
rising  and  falling  water  surfaces,  so  that  the  selected  grid 
should  be  large  enough  to  contain  the  largest  anticipated  volume. 
The  grid  consists  of  active  and  inactive  cells,  the  former  being 
those  at  which  velocities,  temperatures,  and  constituent  con¬ 
centrations  are  computed  and  the  latter  being  those  at  which  no 
computations  are  done.  Inactive  cells  are  identified  by  inputing 
their  widths  as  zero.  The  grid  must  satisfy  the  following  rules: 

’  it  must  be  at  least  two  active  cells  deep  at  every 
segment; 

*  it  must  be  at  least  three  active  cells  long  at  every 
layer ; 
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steps  in  the  waterbody  bottom  profile  are  permitted  as 
.long  as  each  continuous  layer  produced  is  at  least  three 
cells  long; 

the  grid  must  include  a  layer  of  inactive  cells  at  the 
top  and  bottom  and  a  segment  of  inactive  cells  at  the 
left  and  right,  such  that  those  cells  representing 
the  waterbody  are  surrounded  on  all  sides  by  a  layer 
or  segment  of  inactive  cells;  and 

cell  widths  must  not  decrease  from  bottom  to  top  in 
any  segment 

Active  cells  are  those  that  may  contain  water,  even  though  at 
various  times  because  of  flood  events  and  drawdown  they  do  not. 
The  current  boundaries  of  the  computation  on  the  grid  are  de¬ 
fined  by  the  parameters  IL  and  IMAXU1  in  the  longitudinal  and 
KT  and  KMAXM1  in  the  vertical.  The  minimum  values  of  the  para¬ 
meters  IL  and  KT  are  each  two,  indicating  the  waterbody  is  at 
its  maximum  volume.  As  drawdown  occurs,  KT  increases,  so  that 
the  water  surface  is  located  near  layer  KT.  If  necessary,  LARM2 
also  decreases  IL  so  that  segment  I=IL  is  always  at  least  two 
cells  deep. 

Layers  are  identified  by  the  layer  number  K,  with  a  range 
of  1  to  KMAX.  Each  of  these  has  an  elevation  associated  with 
it ,  so  that 

EL  -  h* (KMAX  -  KT+1)  -  Z  +  DTM  (8) 

is  the  water  surface  elevation,  where  DTM  is  the  distance  from 
some  reference  datum  to  the  bottom  of  the  grid.  LARM2  computes 
water  surface  elevations  in  terms  of  the  layer  number,  K,  and 
the  local  elevation,  Z,  rather  than  an  actual  elevation.  Lon¬ 
gitudinal  distances  are  given  in  terms  of  the  segment  number  I. 
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Mapping  these  segments  and  layers  back  onto  the  plan  and  ele¬ 
vation  of  the  reservoir  will  help  to  relate  LARM2  results  to 
specific  reservoir  locations. 

This  mapping  will  also  be  useful  in  establishing  the 
locations  of  inflows  and  outflows.  LARM2  uses  the  following 
definitions: 

upstream  inflows  occur  only  at  I=IL  (left  boundary)  and 
are  defined  by  the  variables  QIN,  TIN,  and  CIN; 

downstream  outflows  occur  only  at  I=IMAXM1  (right  boundary) 
and  are  defined  by  the  variable  QOUT; 

there  must  always  be  an  inflow  and  at  least  one  outflow 
specified,  although  the  associated  flows  may  be  zero; 

tributaries  can  occur  in  any  segment  from  I=IL  to 
I=IMAXM1;  each  tributary  must  have  an  associated  flow, 
temperature,  and  constituent  concentration;  and 

withdrawals  can  be  located  in  any  segment  and  layer  and 
each  must  have  an  associated  flow. 

The  reservoir  orientation  may  be  specified  if  wind  speed 
and  direction  data  are  available  and  if  wind  stress  computations 
are  iesired.  The  general  reservoir  direction  is  determined  by 
the  angle  the  center  line  (x-axis,  positive  to  the  right  down¬ 
stream)  makes  with  north. 

LARM2  as  presently  coded  will  handle  grid  sizes  up  to 
30  segments  long  and  50  layers  deep.  Larger  sizes  will  require 
dimension  statements  to  be  changed.  These  statements  are 
marked  "DIMENS"  in  columns  73  to  80, and  instructions  for  changing 
them  are  located  in  comment  statements  preceding  each  dimension 
statement.  Dimension  statements  are  found  in  the  main  program 
and  subroutines  TRIDAG,  GRIDG,  GRID,  and  GRIDC. 
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4.3  Time- Varying  Boundary  Condition  Data 

In  order  to  model  the  water  and  heat  budgets  accurately 
through  the  simulation  period,  it  is  necessary  to  consider  both 
heat  exchange  at  the  water  surface  and  mass  and  heat  advected  into 
or  out  of  the  waterbody.  These  parameters  are  provided  to  the 
main  code  through  the  subroutine  TVDS,  which  is  an  acronym  for 
time-varying  data  selector.  This  subroutine  performs  two  tasks. 

It  is  called  once  in  order  to  read  the  time-varying  data.  On 
each  succeeding  call,  TVDS  receives  from  the  main  code  the  cur¬ 
rent  simulation  time  and  returns  to  the  main  code  the  value  of 
each  of  the  parameters  specified  in  its  call.  TVDS  must  be  mod¬ 
ified  by  the  user  to  accommodate  his  particular  data.  Minimum 
data  requirements  are  flows  and  temperatures  for  the  upstream 
inflow  (QIN,  TIN)  and  each  tributary  (QTRIB,  TTRIB);  outlet 
flows  at  the  dam  (QOUT)  and  each  withdrawal  (QWD) ;  and  the  heat 
exchange  parameters,  CSHE,  SRO,  and  ET.  If  evaporation  is  to  be 
computed,  the  dew  point  temperature  (TD)  is  required;  wind  speed 
(WA)  and  direction  (PHI)  are  needed  for  wind  stress  computations. 
For  each  constituent  being  simulated,  inflow  and  tributary  con¬ 
centrations  are  also  required.  These  parameters  can  be  supplied 
at  any  time  interval.  Daily  or  more  frequent  data  is  recommended, 
although  in  some  cases  simulations  made  with  constant  parameters 
are  useful . 

The  subroutine  TVDS  is  given  the  simulation  time  (variable 
ELTM)  which  is  measured  in  Julian  days  and  fractions  thereof. 


ELTM  is  computed  in  LARM2  as  TMSTRT  +  N*DLT/86400,  where  N  is 
the  current  iteration  number.  Note  that  an  ELTM  of  1.75  means 
January  1,  6  p.m.  and  that  an  ELTM  of  72.0  is  midnight  March  12/ 
March  13.  A  Julian  date  calendar  is  given  in  Appendix  I. 

In  the  most  general  case,  TVDS  has  available  from  its  orig¬ 
inal  call  several  arrays,  each  one  of  which  contains  time-varying 
data  from  a  separate  disc  or  tape  file.  In  the  example  of 
Dillon  Reservoir,  one  array  contains  all  the  flow  data,  a  second 
array  the  ammonia  concentration  data,  and  a  third  the  pre- 
processed  meteorological  data,  including  an  inflow  temperature 
estimate  for  the  upstream  inflow  and  two  tributaries.  The 
first  two  arrays  contain  daily  data,  and  the  last  one  contains 
data  on  an  irregular  interval.  All  these  files  are  chronolog¬ 
ical,  which  is  essential  for  the  TVDS  algorithm.  Pairs  of  times 
associated  with  each  data  file  are  examined  sequentially  and 
compared  to  ELTM,  beginning  with  the  previously  located  line  of 
data.  When  the  correct  line  is  located,  the  data  is  extracted 
from  the  array  and  set  equal  to  the  proper  subroutine  arguments. 
In  the  example  shown,  all  three  data  arrays  are  examined  in 
this  manner  before  the  time-varying  boundary  condition  data  are 
returned  to  the  main  code.  The  subroutine  also  contains  error¬ 
checking  routines. 

It  is  important  that  the  prospective  user  study  the  TVDS 
subroutine.  Differences  in  data  availability  and  format  make  it 
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easier  for  the  user  to  tailor  TVDS  to  his  data  than  to  attempt 
to  create  a  generalized  TVDS  able  to  handle  all  cases. 


4.4  Initial  Conditions,  Hydraulic  Parameters,  and  Constituent 
Reactions 


Having  established  a  grid  to  represent  the  reservoir  and 
assembled  time-varying  data  for  the  simulation  period,  it  is 
then  necessary  to  select  initial  conditions  to  start  each 
simulation.  These  consist  of  a  water  surface  elevation  and 
temperature  corresponding  to  the  simulation  starting  time. 

This  information  is  coded  on  the  IC  input  card  (Appendix  A). 

If  the  water  quality  constituent  computation  option  is  used, 
initial  conditions  for  each  constituent  are  required.  These 
are  specified  on  the  CC  card. 

The  initial  temperature  (and  constituent)  field  is  speci¬ 
fied  as  a  single  number,  so  that  the  code  begins  with  an  iso¬ 
thermal  field.  With  some  reprogramming,  a  stratified  field  can 
be  used.  For  the  usual  simulation  situation  of  early  spring 
start-up,  the  isothermal  approach  works  well.  When  there 
is  no  available  data  to  begin  a  simulation,  this  may  be  the 
only  possible  approach.  In  flow -dominated  cases,  the  compu¬ 
tations  require  an  initialization  period  approximately  equal  to 
the  reservoir  residence  time  to  damp  out  this  initial  temperature 
estimate.  Therefore,  the  initial  results  up  to  approximately 
a  residence  time  into  the  simulation  may  be  spurious.  If  this 
time  period  is  important,  it  is  suggested  that  the  time-varying 
data  be  duplicated  at  the  start  of  the  simulation  such  that, 
for  example,  March's  meteorological  and  hydrologic  data  is 
followed  by  March  (again),  April,  May,  etc.,  for  a  waterbody 
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with  a  residence  time  of  about  a  month.  This  procedure  may  be 
continued  as  additional  simulations  are  made,  or  results  for 
the  end  of  March  may  be  used  to  define  a  new  initial  elevation 
and  temperature.  This  technique  can  also  be  used  to  extend  a 
previous  simulation  by  selecting  initial  conditions  from  the 
previous  simulation  results. 

The  flow  field  is  always  initialized  to  zero  by  the  code. 
The  solution  technique  is  such  that  the  applied  boundary  con¬ 
ditions  very  quickly  force  an  internally  compatible  flow  field. 

Hydraulic  parameters  consist  of  Ax,  Az>  and  the  Chezy 
coefficient.  These  parameters  are  defined  in  Chapter  2  of 
this  guide,  and  suggested  values  are  coded  in  DATA  statements. 

A  and  A  may  be  varied  by  several  orders  of  magnitude,  with 
the  following  exception.  A  .  the  vertical  dispersion  of  mo- 

mentum  coefficient,  has  limits  of  the  molecular  value 
-6  2 

( 1 . 5X10  m  /s)  as  a  minimum  and  a  value  related  to  the  hori¬ 
zontal  layer  thickness  and  the  time  step  as  a  maximum 
(A  <h2/2At).  The  latter  limit  is  derived  from  the  computation 

Z 

time  step  limit,  Equation  7.  The  value  of  A  varies  spatially 
and  temporally  within  these  limits  and  is  computed  by  the  code 
from  the  Richardson  number,  which  is  the  local  ratio  of  buoy¬ 
ancy  to  shear.  The  value  of  Az  used  in  the  computation  is  a 
function  of  a  base  value  for  this  parameter,  Azo ,  and  the 


Richardson  number: 
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Az  -  Az()  (1  +  10R±)"1/2  (9) 

where  Ri  is  the  Richardson  number.  It  is  necessary  to  set  only 
,  as  the  code  automatically  keeps  the  value  of  within  the 
molecular  and  computational  limits.  Ax  is  presently  set  at  a 
single  value  that  is  constant  over  time  and  space,  but  the 
numerical  formulation  allows  it  to  vary  spatially  and  temporally 
with  additional  coding. 

The  Chezy  coefficient  is  also  space-  and  time-invariant, 
but  can  be  made  variable  with  some  reprogramming.  Boundary 
friction  varies  spatially  with  the  amount  of  side  and  bottom 
area  exposed,  as  well  as  temporally  with  the  flow  field.  Chezy 
may  be  varied  by  a  factor  of  four  from  the  value  given  in  the 
LARM2  DATA  statement . 

The  temperature  equation  dispersion  parameters  may  be 
varied.  These  parameters  are  Dx  and  Dz>  the  longitudinal  and 
vertical  dispersion  coefficients.  Like  A  ,  D  is  temporally 

X  X 

and  spatially  invariant ,  and  its  value  is  controlled  by  the 

FORTRAN  variable  DXI.  Dz  varies  over  time  and  space  and  is 

also  computed  from  the  Richardson  number  and  a  base  value,  Dzo • 

Like  the  momentum  equation  parameter  A  .  D  operates  over  a 

z  z 

range  from  a  minimum  molecular  value  to  a  computational  maximum 
of  h2/2At. 

Finally,  the  distribution  of  solar  radiation  in  the  ver¬ 
tical  is  controlled  by  the  parameters  B  and  y  (FORTRAN  BETA  and 
GAMMA)  as  follows: 
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H  <z)  -  (1  -  8)  H  (z-0)“YZ  (10) 

8  S 

where 

H  solar  radiation  rate  at  depth  z  (°C*m3*s  1 /m2) 

s 

6  fraction  of  H  absorbed  at  the  water  surface 

s 

y  attenuation  rate 

Note  that  these  parameters  control  only  the  distribution  of 
solar  radiation  in  the  vertical  and  that  overall  surface  heat 
exchange,  including  solar  radiation,  is  summarized  in  the  para¬ 
meters  CSHE  and  ET.  If  8=1,  all  the  solar  radiation  is  absorbed 
in  the  current  top  layer  (K=KT)  and  the  value  of  y  is  arbitrary. 

Constituent  Reactions 

For  each  constituent  specified,  constituent  reaction  rates 
and  internal  sources  and  sinks  need  to  be  coded.  For  the  user 
who  desires  only  hydrodynamic  and  temperature  simulations,  this 
section  of  the  user  guide  may  be  skipped.  For  a  detailed  de¬ 
scription  of  the  constituent  transport  computations,  the  user  is 
referred  to  "Developments  in  LARM2 :  A  Longitudinal  —  Vertical ,  Time 
Varying  Hydrodynamic  Reservoir  Model  "  (Edinger  and  Buchak,  1983). 

The  variable  RR(JC,M)  is  the  rate  at  which  constituent  JC 
is  transferred  to  constituent  M  and  may  be  dependent  on  a  number 
of  spatially  and  temporally  varying  variables,  such  as  temperature. 
For  the  case  of  M=JC,  RR  is  negative  and  is  the  linear  decay  rate. 
RR  has  the  units  of  s~ 1 . 

The  variable  HN(I,K)  is  the  spatially  and  temporally  varying 
source  and  sink  term  for  constituent  JC.  HN(I,K)  has  the  units 
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mg  1  •  m  • s  .  The  user  is  required  to  augment  this  array  by 

the  amount  of  material  added  to  or  subtracted  from  constituent 
JC  due  to  the  reactions  specified  by  RR. 

The  user  must  examine  DO  loop  1840  in  the  LARM2  code  and 
modify  it  to  compute  RR  and  HN  as  required  for  his  constituents. 
The  user  may  also  specify  different  sources  and  sinks  in  the 
surface  or  bottom  layers  by  separate  statements  for  layers 
K=KT  or  K=KB.  The  computation  for  each  constituent  concentration 
(DO  loop  1820)  begins  with  an  initialization  of  the  HN  array. 
Then,  for  each  location  I,K,  the  reaction  rate  (RR)  array  is  com¬ 
puted,  one  FORTRAN  statement  for  each  constituent  pair  JC,M. 

HN  is  then  augmented  by  the  rate  at  which  constituent  JC  is 
transferred  to  or  from  constituent  M.  The  remainder  of  the 
computation  requires  no  user  c ianges  and  includes  the  augmenting 
of  HN  by  external  sources  and  sinks,  the  retrieving  of  the  trans¬ 
port  coefficient  vectors,  and  the  call  to  the  subroutine  TRIDAG 
for  layer-by-layer  solution  of  the  transport  equation. 
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4.5  Input,  Output  and  Computer  Resource  Requirements 

The  input  data  stream  is  described  in  Appendix  A,  and  an 
example  data  set  for  Dillon  Reservoir  is  given  in  Appendix  B. 
LARM2  output  for  this  data  set  is  given  in  Appendix  C.  The 
output  shown  there  is  for  IFORM=0,  that  is,  compressed  (for 
80-character-width  printers)  and  shortened  results.  The  normal 
form  of  the  output  (for  132-character-width  printers)  can  be 
obtained  by  setting  the  parameter  IF0RM=1.  Although  an  example 
of  this  more  complete  output  is  not  given,  the  user  may  quickly 
make  his  own  by  changing  the  IFORM  parameter  (PR  card,  field  1) 
from  0  to  1.  The  additional  information  generated  by  the  longer 
output  is  as  follows: 

‘  the  long  form  produces  three  geometry  tables:  a  grid 
showing  active  and  inactive  cells,  a  grid  showing  cell 
widths  as  augmented  by  LARM2  (widths  assigned  to  in¬ 
active  cells),  and  an  elevation-area-volume  table; 

’  the  initial  conditions  for  Z2,  U,  YF,  and  T2  (and  C2): 
the  velocity  fields  are  always  initialized  to  zero  at 
the  start  of  the  computation,  Z2  and  T2  are  initialized 
through  the  IC  card; 

eight  additional  segment  results  for  Z2,  V,  W,  and  T2 
(and  C2),  using  the  additional  columns  available  on 
11-inch  by  14-inch  paper. 

A  second  type  of  output  is  available  from  the  LARM2  code, 
and  it  is  useful  in  postprocessing  applications  such  as  the 
vector  plotting  code  W.  Cell  coordinate  positions  are  written 
unformatted  to  TAPE61,  followed  by  U  and  IF  fields  at  times 
selected  on  the  PL  card  (see  Appendix  A).  This  information  is 
used  by  W  to  plot  displacement  vectors  on  a  scaled  grid. 


However,  by  saving  additional  information  on  TAPE61,  time  series 
plots  of  other  variables  may  be  obtained,  as  well  as  statistical 
analyses  of,  for  example,  outflow  temperatures  or  velocity 
records  at  a  particular  reservoir  location. 

The  code  and  test  data  as  supplied  requires  254,000  octal 
words  of  storage  on  a  Controlled  Data  Corporation  (CDC)  175. 

This  storage  requirement  depends  on  both  the  grid  size  and 
the  amount  of  time-varying  boundary  condition  data  stored. 

LARM2  requires  0.00026  s  of  central  processor  (cpu)  time  for 
each  time  step  iteration  for  each  active  cell  at  a  high  level 
of  compiler  optimization.  The  Dillon  Reservoir  example  uses 
90.84  s  of  cpu  time  for  its  366  active  cells  for  a  simulation 
of  960  steps,  plus  9.5  s  for  compilation. 
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4.6  Dillon  Reservoir  Example 

Dillon  Reservoir  is  located  in  Colorado  approximately  90  km 
west  of  Denver  in  the  Rocky  Mountains.  The  reservoir  is  used  for 
water  supply  as  well  as  for  receiving  effluent  from  a  wastewater 
reclamation  facility  serving  nearby  ski  resorts.  The  problem 
addressed  with  the  LARM2  simulations  is  the  dilution  and  decay 
of  ammonia  from  the  treatment  plant  located  on  a  tributary  and 
from  nonpoint  sources  in  the  basin.  The  vertical  and  longitudinal 
distribution  of  ammonia  in  the  reservoir  arms  is  of  particular 
interest  since  several  species  of  fish  migrate  through  the  arms 
to  the  shallow  tributaries  to  spawn.  The  unusual  geometry  and 
the  inclusion  of  constituent  computation  make  this  a  complex 
case.  However,  the  user  will  find  elements  of  many  applications 
in  the  Dillon  Reservoir  example  and  is  encouraged  to  study  it. 

Figure  3  shows  Dillon  Reservoir  in  plan.  In  order  to  in¬ 
clude  both  the  Blue  River  and  Tenmile  Creek  arms  in  the  LARM2  grid, 
the  center  line  was  drawn  as  shown  in  Figure  4.  Cross-sectional 
data  at  500-ft  intervals  were  obtained  and  formatted  for  GEDA. 

The  results  of  the  GEDA  run  were  satisfactory  with  regard  to  the 
comparison  of  the  design  report  elevation-area- volume  table,  so 
GIN  was  run  with  the  GEDA  results  as  input.  Several  GIN  runs 
were  required  to  produce  an  acceptable  geometry;  some  modi¬ 
fication  of  the  widths  was  necessary  to  match  the  top  layer 
surface  area  with  the  design  report  surface  area  ar»d  produce  a 
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good  volume  fit  elsewhere  in  the  vertical.  The  final  GIN  run 
provided  the  BA  cards  for  use  as  LARM2  input  (Appendix  B). 

The  final  grid  is  shown  in  Figures  2  and  4.  The  output  from 
GIN  also  gives  the  information  required  on  the  GE  input  card, 
namely  the  grid  size,  Ax  and  h,  and  the  elevation  of  the  grid 
relative  to  the  reference  datum. 

Coincident  with  establishing  a  grid  representation  of  Dillon 
Reservoir,  inflows  and  withdrawals  and  inflow  temperatures  and 
ammonia  concentrations  were  identified.  For  the  grid  shown  in 
Figure  2,  the  following  nomenclature  was  adopted: 

Blue  River:  QIN,  TIN,  CIN(l)  -  upstream  inflow 

Snake  River:  QTRIB(l) ,  TTRIB(l),  CTRIB(1,1)  -  segment  9 

Tenmile  Creek:  QTRIB(2),  TTRIB<2),  CTRIB(2,1)  -  segment  17 
Roberts  Tunnel:  QWD(l)  -  segment  9,  layer  28 
Dam:  QWD(2)  -  segment  10,  layer  35 

Note  that  because  of  the  bending  of  the  center  line  up  the  Ten- 
mile  Creek  arm,  the  dam  outflow  is  now  taken  as  a  withdrawal, 
not  an  outlet  as  defined  for  LARM2.  The  code  requires  at  least 
one  outlet,  so  one  is  specified  (on  the  OU  card) ,  but  its  flow  is 
always  set  to  zero  in  TVDS  and  its  layer  specification  on  the 
OU  card  is  arbitrary.  Also  note  that  the  inflow  as  defined  for 
LARM2  is  the  Blue  River  and  that  Tenmile  Creek  is  treated  as  a 
tributary,  that  is,  it  is  not  considered  a  momentum  source  for 
the  reservoir.  The  assignment  of  these  inflows  and  withdrawals 
allows  completion  of  the  TR,  TL,  WD,  WI,  and  WK  cards. 
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The  time-varying  data  were  made  available  in  three  files  for 
flows,  inflow  ammonia  concentrations,  and  meteorological  data. 

The  meteorological  data  were  preprocessed  from  a  record  con¬ 
sisting  of  date  and  time,  air  and  dew  point  temperatures,  wind 
speed  and  direction,  and  observed  cloud  cover  to  a  record  consist¬ 
ing  of  year,  Julian  date,  dew  point  temperature,  wind  speed  and 
direction,  equilibrium  temperature,  coefficient  of  surface  heat 
exchange,  estimated  solar  radiation,  and  response  temperature. 
Response  temperature  is  computed  as 


k  (Tr 


E) 


where 


(ID 


assumed  depth  (m) 
response  temperature  (°C) 
time  (s) 

kinematic  coefficient  of  surface  heat 
exchange  (m/s) 

equilibrium  temperature  (°C) 

The  response  temperature  was  used  to  estimate  TIN,  TTRIB(l),  and 
TTRIB(2)  because  inflow  temperature  records  were  not  available. 
Equation  11  represents  a  heat  balance  for  a  fully-mixed  condi¬ 
tion  in  which  change  in  heat  storage  is  balanced  by  surface 
heat  exchange.  By  adjusting  the  depth,  D,  to  fit  the  few  avail¬ 
able  observations,  an  estimate  of  the  seasonal  inflow  tempera¬ 
ture  was  made. 


r 

t 

k 

E 


The  subroutine  TVDS  was  constructed  around  the  three 


time-varying  data  files.  On  LARM2's  first  call  to  TVDS  each 
of  the  three  files  is  read,  converted  to  SI  units,  and  stored 
in  arrays.  On  every  other  call  to  TVDS,  the  arrays  are  searched 
to  locate  the  time- varying  data  corresponding  to  the  LARM2 
time  variable  ELTM.  Since  the  data  files  are  sequential,  the 
lower- loop  parameter  can  be  updated  every  time  the  correct  data 
is  located.  This  updating  permits  the  next  search  to  begin  at 
the  previous  successful  location. 

Observed  data  for  the  1975  study  year  consisted  of  daily 
water  surface  elevations  and  weekly  temperature  profiles  at  a 
single  location.  Initial  conditions  for  each  simulation  were 
selected  from  these  data.  The  first  simulations  were  for  a  few 
days  only  and  were  checked  for  correct  time-varying  data  and 
reasonableness  of  results.  Long-term  simulations  were  then 
made  beginning  with  January  2  and  continuing  throughout  the 
year.  It  was  found  that  a  At  of  900  s  was  too  large  to  accommo¬ 
date  the  velocities  of  the  fall  turnover,  so  that  a  At  of  600  s 
was  used  thereafter  for  each  long-term  simulation.  The  results 
of  long-term  simulations  were  compared  to  the  weekly  temperature 
observations.  These  comparisons  showed  the  LARM2  computed  tem¬ 
peratures  too  high  in  the  spring  and  too  low  in  the  fall.  A 
heat  budget  analysis  showed  that  the  error  was  due  to  inflow 
temperatures  that  responded  too  quickly  to  the  annual  heating¬ 
cooling  cycle  due  to  meteorological  conditions.  By  increasing 
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the  parameter  D  in  Equation  11,  a  better  inflow  temperature 
estimate  was  obtained.  The  LARM2  simulations  using  this  new 
record  produced  satisfactory  agreement  between  computed  and 
observed  temperatures  within  the  limits  of  extrapolating  me¬ 
teorological  data  from  a  distant  site. 

Weekly  velocity  plots  were  also  examined  using  the  post¬ 
processor  W.  Although  no  current  measurements  were  available 
for  Dillon  Reservoir,  the  plots  were  studied  for  overflow,  inter¬ 
flow,  and  underflow  in  the  reservoir  arms  as  well  as  for  be¬ 
havior  during  turnover. 

Once  the  computed  velocity  and  temperature  fields  were 
shown  to  be  realistic,  ammonia  simulations  were  begun.  Once 
again,  data  from  the  observation  record  were  used  for  initial 
conditions  and  verification.  The  first  simulations  proved  to 
be  satisfactory.  Finally,  new  inflow  ammonia  concentrations, 
reflecting  an  expected  increase  in  treatment  plant  loading, 
were  substituted  for  the  existing  ammonia  concentration  file. 

The  results  of  this  simulation  were  examined  for  effects  on 
the  fish  community  with  regard  to  potentially  toxic  levels  and 
possible  interference  with  spawning  activities. 

The  time-varying  data  files  given  in  the  test  data  set  and 
shown  in  Appendix  B  are  for  two  weeks  of  the  1975  complete  data 
set.  The  example  simulation  shown  in  Appendix  C  used  results 
from  the  full-year  simulation  to  obtain  initial  conditions. 

For  the  two  weeks  of  data  shown,  a  At  of  900  s  proved  to  be 
acceptable. 


REFERENCES 


Edinger,  J.E.,  Brady,  D.K.,  and  Geyer,  J.C.  (1974),  Heat  Exchange 
and  Transport  in  the  Environment ,  Publication  No.  74-049-00-3, 
Electric  Power  Research  Institute,  Palo  Alto,  Calif. 

Edinger,  John  Eric  and  Edward  M.  Buchak  (1979),  A  Hydrodynamic, 
Two-Dimensional  Reservoir  Model:  Development  and  Test  Appli¬ 
cation  to  Sutton  Reservoir  Elk  River,  West  Virginia,  Contract 
No.  DACW27-76-C-0089  ,  U.  S.  Army  Engineer  Division,  Ohio  River 
Cincinnati. 

Edinger,  John  Eric  and  Edward  M.  Buchak  (1983),  Developments  in 
LARM2 :  A  Longitudinal-Vertical,  Time-Varying  Hydrodynamic 

Reservoir  Model,  Technical  Report  in  press,  prepared  by  J.  E. 
Edinger  Associates  for  U.  S.  Army  Engineer  Waterways  Experi¬ 
ment  Station,  Vicksburg,  Miss. 

Hydrologic  Engineering  Center  (1976),  Geometric  Elements  from 
Cross  Section  Coordinates.  Davis,  Calif. 


QTPIB 


FIGURE  3  Dillon  Reservoir 


FIGURE  4  Dillon  Reservoir  LARU  Segmentation 


GENERAL 


Descriptions  for  each  of  the  cards  (or  card  images)  used  as 
input  to  LARM2  are  presented  on  the  following  pages.  There 
are  three  data  types: 

alphanumeric; 


real  (may  be  left  or  right  justified);  and, 
integer  (must  be  right  justified). 

There  are  eleven  data  fields  per  card : 


Field 

0 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 


Length 

2 

6 

8 

8 

8 

8 

8 

8 

8 

8 

8 


Columns 

1-2 

3-8 

9-16 

17-24 

25-32 

33-40 

41-48 

49-56 

57-64 

65-72 

73-80 


This  format  is  similar  to  that  used  in  the  Hydrologic  Engineering 
Center  codes. 
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TITLE  CARDS  (REQUIRED) 

Title  cards  may  be  used  to  identify  applications,  simulations  or 
other  parameters  as  required.  These  cards  do  not  otherwise  affect 
the  computation. 


CARD  T1 


The  text  input  with  this  card  appears  in  the  header  with  each  ve¬ 
locity,  temperature,  and  constituent  field  snapshot,  as  well  as  in 
the  overall  header  printed  once  for  each  simulation. 


Field 

Parameter 

Value 

Description 

0 

AID 

T1 

alphanumeric: 

card  identification 

1-9 

TITLE(l.J), 
J=1 , 7 

alphanumeric : 
fication  text 

simulation  identi- 

CARD  T2 


The  text  input  using  this  card  always  appears  in  the  overall  simu¬ 
lation  header  and  appears  in  the  first  snapshot  header  if  IFORM=l . 


Field 

Parameter 

Value 

Description 

0 

AID 

T2 

alphanumeric:  card  identification 

1-9 

CARD  T3 

TITLE(2, J), 
J=1 , 7 

alphanumeric:  simulation  identi¬ 
fication  text 

The  text  input  using  this  card  always  appears  in  the  overall  simu¬ 
lation  header  and  appears  in  the  first  snapshot  header  if  IF0RM=1. 


Field 

Parameter 

Value 

Description 

0 

AID 

T3 

alphanumeric: 

card  identification 

1-9 

TITLE( 3 , J) , 
J=1 , 7 

alphanumeric : 
fication  text 

simulation  identi- 

CARD  T4 


The  text  input  using  this  card  always  appears  in  the  overall  simu¬ 
lation  header  and  appears  in  the  first  snapshot  header  if  IF0RM=1. 


Field 

Parameter 

Value 

Description 

0 

AID 

T4 

alphanumeric: 

card  identif icatioi 

1-9 

TITLE(4, J) , 
J=1 , 7 

alphanumeric: 
fication  text 

simulation  identi- 

Note:  The  variable  TITLE  is  dimensioned  to  accommodate  seven  ten- 
character  words  per  line,  which  is  characteristic  of  CDC 
processors.  The  number  of  characters  per  word  varies:  IBM 
accommodates  four,  Univac  six.  In  adapting  LARM2  to  other 
machines,  both  the  dimension  of  TITLE  and  FORMAT  statements 
500,  610,  615,  764,  765,  and  9000  may  need  to  be  changed. 
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GEOMETRY  CARD  (REQUIRED) 
CARD  GE 


The  geometry  card  specifies  the  grid  size  (number  of  segments  and 


layers),  the  length  and  height  of 
location  of  the  grid  relative  to 


Field 

Parameter 

Value 

0 

AID 

GE 

1 

IMAX 

>5 

2 

KMAX 

>4 

3 

DLX 

+ 

4 

HIN( 1,1) 

+ 

5 

DTM 

-,o,+ 

the  individual  cells,  and  the 
a  reference  datum. 

Description 

alphanumeric:  card  identification 

integer:  number  of  segments 

integer:  number  of  layers 

real:  Ax,  cell  length,  m 

real:  h,  cell  height,  m 

real:  distance  from  reference 

datum  to  lower  grid  boundary,  m 


Note:  LARM2  requires  a  grid  with  IMAX>5  and  KMAX>4.  A  typical 

grid  has  IMAX  in  the  10-to-30  range,  KMAX  in  the  15-to-50 
range.  Grid  sizes  larger  than  30*50  require  changing  the 
code's  dimension  statements.  Instructions  for  doing  so  are 
given  in  comment  cards  preceding  each  dimension  statement , 
each  of  which  is  marked  "DIMENS"  in  columns  73-80. 

Typical  cell  sizes  are  DLX  in  the  0. 5-km-to-10-km  range  and 
HIN(1 , 1)  in  the  0.5-m-to-3-m  range. 


TIME  CARD  (REQUIRED) 


CARD  TM 

This  card  determines  the  simulation  length  by  specifying  the  tem¬ 
poral  iteration  size  and  the  number  of  iterations. 


Field 

Parameter 

Value 

Description 

0 

AID 

TM 

alphanumeric:  card  identification 

1 

DLT 

+ 

real:  At  computation  time  step,  s 

2 

NSTEPS 

+ 

integer:  number  of  iterations 

Note:  The  simulation  length  is  NSTEPS  *  DLT,  in  seconds. 
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PLOT  CARD  (REQUIRED) 
CARD  PL 


This  card  specifies  whether,  and  at  what  frequency,  velocity  fields 
are  saved  for  postprocessing  by  W,  the  vector  plotting  code. 


Field  Parameter  Value 

0  AID  PL 

1  IPLOT  0 

1 

2  Ml  0,  + 

3  M2  0,+ 


Description 

alphanumeric:  card  identification 

integer:  no  information  saved 

integer:  information  saved 

integer:  information  is  saved  be¬ 

ginning  at  iteration  Ml;  if  IPLOT=0, 
this  paramter  has  no  effect 

integer:  information  is  saved 

every  M2  iterations;  if  IPL0T=0, 
this  parameter  has  no  effect 


Note:  If  IPLOT=l ,  Ml=960,  and  M2*96  for  example,  postprocessing  in¬ 

formation  is  saved  every  96  iterations,  beginning  at  itera¬ 
tion  960. 
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PRINT  CARD  (REQUIRED) 
CARD  PR 


This  card  specifies  the  form  of  the  printed  output,  as  well  as  the 
frequency  of  the  velocity,  temperature,  and  constituent  field 
snapshots. 


Field 

Parameter 

Value 

Description 

0 

AID 

PR 

alphanumeric:  card  identification 

1 

I  FORM 

0 

integer:  output  shortened  (cell 
widths  and  volurae-area-elevation 
table  suppressed)  and  compressed 
to  fit  an  8J-inch  page  width 

1 

integer:  normal  output  using  en¬ 

tire  132-column  page  width  (14 
inches) 

2 

N1 

o,+ 

integer:  snapshots  are  printed 

beginning  at  iteration  N1 

3 

N2 

o,+ 

integer:  snapshots  are  printed 
every  N2  iterations 

The  parameters  N1  and  N2  are  used  like  the,  parameters  Ml  and 
M2  on  the  PL  card.  The  segments  for  which  elevations,  veloc¬ 
ities,  and  temperatures  are  printed  in  the  snapshots  are 
specified  on  the  10,  II,  and  12  cards. 


Note: 
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SEGMENT  PRINT  SELECTOR  CARDS  (REQUIRED) 


These  cards  permit  the  user  to  select  which  segment  results  are  to 
be  printed  for  IFORM=0  and  IF0RM=1.  Nine  segments  may  be  selected 
for  IFORM=0  (short  form)  and  17  segments  for  IF0RM=1  (long  form). 


CARD  10 

Field 

Parameter 

Value 

Description 

0 

AID 

10 

alphanumeric : 

card  identification 

1-9 

10(1), 

1-1,9 

1< I0< IMAX 

integer:  segment  number 

CARD  11 

Field 

Parameter 

Value 

Description 

0 

AID 

11 

alphanumeric : 

card  identification 

1-10 

11(1), 

1*1,10 

1< Il< IMAX 

integer:  segment  number 

CARD  12 

Field 

Parameter 

Value 

Description 

0 

AID 

12 

alphanumeric: 

card  identification 

1-7 

IKD, 

1*11,17 

1<I1<IMAX 

integer:  segment  number 

Note:  If 

IMAX<9  or 

IMAX<17 ,  segment  numbers  may 

be  repeated  to  com- 

plete  the  card  requirements. 


INITIAL  CONDITION  CARD  (REQUIRED 


N 


CARD  IC 


This  card  defines  the  starting  time  in  Julian  days  and  the  corre¬ 
sponding  initial  water  surface  elevation  and  temperature. 


Field 

Parameter  Value 

Description 

0 

AID 

IC 

alphanumeric:  card  identification 

1 

TMSTRT 

0<TMSTRT<366 

real:  simulation  starting  time, 

Julian  days 

2 

KT 

KMAX-2<KT<2 

integer:  layer  number  for  initial 

water  surface  elevation 

3 

ZI 

-0. 8*HIN( 1 , 1 )<ZI 
<0 . 25*HIN( 1,1) 

real:  initial  water  surface  ele¬ 

vation  relative  to  layer  KT,  m 
(positive  down) 

4 

TI 

0<TI<100 

real:  initial  temperature,  °C 

Note: 

KT  and  ZI 

give  the  water 

surface  elevation  and  TI  the  tempera 

ture  corresponding  to  the  time  TMSTRT.  Values  of  3  for  KT 
and  0.1  m  for  ZI,  for  example,  place  the  initial  water  sur¬ 
face  elevation  0.1  m  below  the  top  of  layer  3.  The  initial 
temperature  field  is  isothermal  at  TI.  The  simulation  begins 
at  TMSTRT  and  ends  at  TMSTRT  +  NSTEPS  *  DLT/86400,  in  days. 


All 


CONSTITUENT  DEFINITION  CARD  (REQUIRED) 
CARD  CD 


This  card  specifies  the  number  of  water  quality  constituents  to  be 


simulated  and  also  permits  these 
off  for  a  particular  simulation. 

Field  Parameter  Value 

0  AID  CD 


transport  computations  to  be  turned 
Description 

alphanumeric:  card  identification 


1  NC  0<NC<10  integer:  number  of  constituents 

to  be  simulated;  if  NC>0,  the  fol¬ 
lowing  card,  CONSTITUENT  INITIAL 
CONCENTRATION  CARD,  is  required 


2  ICC 


0 


integer:  constituent  transport 

computations  suppressed 


1 


integer:  constituent  transport 

computations  performed 


Note:  The  user  generally  will  complete  a  hydrodynamic  and  tempera¬ 
ture  simulation  to  his  satisfaction  (ICC=0)  before  going  on 
to  the  water  quality  constituent  computations  (ICC=1). 

Water  quality  reaction-iteration  rates  also  must  be  specified 
for  each  water  quality  constituent.  Since  these  rates  are 
usually  a  function  of  temperature,  they  must  be  computed  in 
LARM2 .  The  location  for  the  user  to  insert  this  computation 
is  discussed  in  Section  4.4  of  this  guide. 


CARD  CC 


This  card  contains  the  initial  concentration  for  each  of  the  NC 
constituents  and  is  required  if  NC>0. 


Field  Parameter  Value 

0  AID  CC 

1-10  CI(JC)  CI(JC)>0 

JC=1 ,NC 


Description 

alphanumeric  card  identification 

real:  initial  concentrations  for 

each  of  NC  constituents 


Note:  LARM2  initializes  the  entire  reservoir  to  CI(JC). 


METEOROLOGICAL  PARAMETER  CARD  (REQUIRED) 


CARD  MP 


This  card  specifies  the  reservoir  orientation  for  wind  stress  com¬ 
putations  and  also  whether  evaporation  is  included  in  the  water 
budget  computations. 

Field  Parameter  Value  Description 

0  AID  MP  alphanumeric:  card  identification 


1  PHIO  0<PHIO<2tt  real:  the  angle  between  the  posi¬ 

tive  x-axis  (reservoir  centerline 
in  the  flow  direction)  and 
north,  rad 

2  IEVAP  0  integer:  evaporation  rates  not 

computed  and  not  included  in  water 
budget  computations 


1  integer:  evaporation  rates  computed 

and  included  in  water  budget  compu¬ 
tations 


Note:  For  a  reservoir  flowing  from  east  to  west,  PHIO=ir/2;  southwest 

to  northeast,  PHI0=5ir/4.  This  convention  is  similar  to  the 
meteorologist's  convention  for  wind  direction. 

Evaporation  rates  are  sometimes  given  in  the  inflow  record  as 
QINnet=QIN-QET.  This  is  the  case  for  which  IEVAP  should  be 
set  to  0.  The  contribution  of  evaporation  to  surface  heat 
exchange  is  always  considered,  regardless  of  the  value  of 
IEVAP. 
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OUTLET  CARD  (REQUIRED 


CARD  OU 


This  card  specifies  the  number  of  reservoir  outlets  and  their  ele¬ 
vations.  At  least  one  outlet  is  required. 

Field  Parameter  Value  Description 

0  AID  OU  alphanumeric:  card  identification 

1  NOUT  1<N0UT<9  integer:  number  of  outlets  lo¬ 

cated  at  right  grid  boundary 

2-10  KOUT(J),  2<K0UT(J)  integer:  layer  number  for  each 

J=1 ,NOUT  <KMAX-1  outlet 


Note:  For  each  outlet  specified,  a  flow  needs  to  be  supplied  by  the 
subroutine  TVDS.  At  least  one  outlet  is  required,  although 
its  flow  may  be  set  to  zero. 

A  single  outlet  spanning  a  number  of  layers  may  be  specified 
by  breaking  the  single  outlet  into  several  and  dividing  the 
outflow  among  them. 
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TRIBUTARY  CARD  (REQUIRED) 
CARD  TR 


This  card  specifies  the  number  of  reservoir  tributaries. 
Field  Parameter  Value  Description 


0 

1 


AID  TR  alphanumeric:  card  identification 

NTRIB  0<NTRIB<10  integer:  number  of  reservoir  trib¬ 

utaries;  if  NTRIB>0,  then  the  fol¬ 
lowing  card,  TRIBUTARY  LOCATION 
CARD,  is  required 


Note:  For  each  tributary,  an  inflow  and  inflow  temperature  need  to 

be  supplied  by  the  subroutine  TVDS. 
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TRIBUTARY  LOCATION  CARD  (OPTIONAL 


CARD  TL 


This  card  gives  the  reservoir  segment  numbers  at  which  tributary 
inflows  enter  and  is  required  if  NTRIB>0. 

Field  Parameter  Value  Description 

0  AID  TL  alphanumeric:  card  identification 

1-10  ITRIB(J),  2<ITRIB( J)  integer:  segment  number  at  which 

J=1,NTRIB  <IMAX-1  each  tributary  enters 


Note:  Tributary  flows  are  entered  in  the  segment  specified  here 

and  at  a  layer  that  corresponds  to  their  density.  Tributaries 
above  the  current  upstream  boundary  segment,  IL,  on  any  par¬ 
ticular  time  step  are  combined  with  the  upstream  inflow. 
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WITHDRAWAL  CARD  (REQUIRED) 
CARD  WD 


This  card  specifies  the  number  of  withdrawals  from  the  reservoir. 

Field  Parameter  Value  Description 

0  AID  WD  alphanumeric:  card  identification 

1  NWD  0<NWD<10  integer:  number  of  reservoir  with¬ 

drawals;  if  NWD>0,  then  the  follow¬ 
ing  cards,  WI  and  WK,  are  required 


Note:  For  each  withdrawal,  an  outflow  rate  needs  to  be  supplied  by 
subroutine  TVDS. 
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WITHDRAWAL  SEGMENT  CARD  (OPTIONAL) 


CARD  WI 


This  card  specifies  the  longitudinal  location  (segment)  for  each 
withdrawal  and  is  required  if  NWD>0. 

Field  Parameter  Value  Description 

0  AID  WI  alphanumeric:  card  identification 

1-10  IWD(J),  2<IWD(J)  integer:  segment  number  for  each 

J=1 , NWD  < IMAX-1  of  NWD  withdrawals 

Note:  The  layer  number  for  each  withdrawal  is  specified  on  the  WK 
card,  which  follows. 
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This  card  specifies  the  vertical  location  (layer)  for  each  with¬ 
drawal  and  is  required  if  NWD>0. 

Field  Parameter  Value  Description 

0  AID  WK  alphanumeric:  card  identification 

1-10  KWD(J),  2<KWD(J)  integer:  layer  number  for  each 

J=1 ,NWD  <KMAX-1  of  NWD  withdrawals 


Note:  The  segment  number  for  each  withdrawal  is  specified  on  the 
preceding  card.  Each  individual  withdrawal  must  have  its 
segment  and  layer  specified  in  the  same  field  on  the  WI  and 
WK  cards,  respectively. 

A  single  withdrawal  spanning  a  number  of  layers  may  be 
specified  by  breaking  it  into  several.  Withdrawals  above 
the  current  upstream  boundary  segment,  IL,  or  above  the 
current  water  surface  layer,  KT,  are  ignored. 
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BATHYMETRY  CARDS  (REQUIRED 


CARDS  BA 


These  cards  contain  the  cell  widths  for  each  of  the  IMAX  KMAX 
cells,  which  are  taken  from  cross-section  drawings  and  represent 
an  average  value  for  each  cell.  GEDA  produces  widths  as  a 
function  of  elevation;  GIN  (see  Appendix  F)  uses  GEDA  output  and 
produces  BA  cards  used  by  LARM2. 

Field  Parameter  Value  Description 

0  AID  BA  alphanumeric:  card  identification 

1-10  B( I ,K)  B(I,K)>0  real:  cell  widths,  m  (left  to 

right,  top  to  bottom,  filling  each 
BA  card) 


Note:  For  the  grid  boundary  segments  (1=1  and  I=IMAX)  and  boundary 

layers  (K=l  and  K=KMAX),  B  must  be  set  to  zero;  inactive  cells 
also  are  identified  by  setting  their  widths  to  zero. 

The  read  statement  for  these  cards  is 
READ(5, 500)  ( (B( I ,K) , K=1 ,KMAX) , 1=1 , IMAX) . 
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These  cards  contain  the  time-varying  boundary  condition  data  for  the 
simulation.  Since  they  are  read  by  the  user-written  subroutine  TVDS, 
no  format  is  specified  here.  A  discussion  of  the  time-varying  bound¬ 
ary  condition  data  and  TVDS  can  be  found  in  Section  4.3. 


T1USER  GUIDE  EXAMPLE  —  OILLON  RESERVOIR,  COLORADO 
T2UPSTREAM  INFLOW,  TWO  TRIBUTARIES,  TWO  WITHDRAWALS 
T30NE  CONSTITUENT  <  AMMONIA) 

T4IFORM=0  FOR  SHORTENED  AND  COMPRESSED  OUTPUT 


—  2/18/81 
,  AND  NO  OUTFLOW 
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62.8  52.9  52.2 

51.9  39.2  52.2 
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53.4  37.5  52.2 

51.0  36.8  52.2 

49.3  28.9  52.2 

95.3  90.3  52.4 


B5 


75. 

224.58 

37.0 

75. 

224.71 

40.0 

75. 

225.09 

43.0 

75. 

225.21 

44.0 

75. 

225.24 

44.0 

75. 

225.46 

42.0 

75. 

225.58 

42.0 

75. 

225.71 

40.0 

75. 

225.84 

42.0 

75. 

225.97 

40.  0 

75. 

226.08 

36.0 

75. 

226.21 

34.0 

75. 

226.33 

34.0 

75. 

226.45 

37.0 

75. 

226.58 

41.0 

75. 

226.71 

37.0 

4.6 

34.0 

3620.1 

5.8 

11.0 

1065.9 

0.0 

0.0 

0.0 

2.3 

27.0 

0.0 

0.0 

0.0 

930.5 

0.0 

0.0 

2754.3 

8.1 

33.0 

1745.4 

11  .5 

31  .0 

2380.5 

1  .2 

36.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

2.3 

27.0 

0.0 

1.2 

36.0 

1665.7 

4.6 

33.0 

4266.7 

3.5 

20.0 

1912.1 

11.5 

33.0 

1540.3 

82.3 

88.5 

52.4 

70.1 

57.2 

52.5 

50.5 

35.7 

52.4 

52.6 

35.7 

52.4 

56.6 

53.3 

52.4 

68.8 

80.3 

52.4 

88.6 

63.1 

52.4 

118.7 

65.4 

52.5 

50.6 

35.4 

52.5 

49.1 

33.2 

52.4 

47.4 

28.2 

52.4 

48.4 

25.7 

52.4 

57.8 

64.9 

52.4 

84.7 

91.8 

52.5 

66.4 

66.5 

52.5 

105.8 

54.7 

52.5 

«  .  ■  .4 


Appendix  C 
Example  Output  Data 


L  A  R  H  —  LATERALLY  AVERAGED  RESERVOIR  MODEL 

VERSION  TWO  ~  MARCH  1981  --  WAYNE*  PENNSYLVANIA 

USER  GUIDE  EXAMPLE  —  DILLON  RESERVOIR.  COLORADO  —  2/18/81 
UPSTREAM  INFLOW*  TWO  TRIBUTARIES*  TWO  WITHDRAWALS*  AND  NO  OUTFLOW 
ONE  CONSTITUENT  (AMMONIA) 

IFORMrO  FOR  SHORTENED  AND  COMPRESSED  OUTPUT 


CONTROL  PARAMETERS 

IMAX=  18  KMAX-  36  DLX=  552.75  M  HIN(1.1>=  2.00  M  DTH=2679.30  M 
DLT=  900.00  S  NSTEPSs  A8 


IPLOTs 

0  Ml  = 

0 

M2  = 

48 

I FORM= 

0  Nl  = 

0 

N2  = 

48 

10  = 

1  2 

3 

4 

5  7 

9  12  15 

11  = 

1  2 

3 

4 

5  6 

7  8  9  10 

11  (12 

CARO  >  = 

11 

12 

13  14 

15  16  17 

TMSTRT 

=  214. 

50  DAYS 

KT=  3 

ZI=  -1.2500  M 

NC= 

1  ICC 

- 

1  WITH  INITIAL  VALUES  OF 

PHI0  = 

3.90 

RAO 

IEVAP=1 

NOUT  = 

1  AT 

K= 

10 

NTRIB= 

2  AT  1  = 

9  17 

NWD= 

2  AT 

1  = 

9 

10 

AND  AT 

K  = 

28 

35 

INTERNALLY-SET  PARAMETERS 

AX=  . 1000E*01  M*M/S 
AZMIN=  .1500E-05  M*M/S 
AZ0=  .3000E-03  M*M/S 
AZMAX=  .1889E-02  M*M/S 
BETAS  1.0000 
CHZY=  70.0000  SORT < H) /S 
0X1=  .1000E+01  M*M/S 
0ZMIN=  .1400E-04  M* M/S 
DZ0=  .2800E-03  M*M/S 
DZMAX=  .1889E-02  M*M/S 
GAMMA=  .4900  1/M 
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LAR42  3/81  JESS  USER  GUIDE  EXAMPLE  —  DILLON  RESERVOIR*  COLORAOO  —  2/18/81 

ELTM=  218.00000  DAYS  <N  =  48)  1975 

BIN-  4.13  CMS  TIN=  9.8  C 

QOUT  ( CMS ) -  0.00 

KOUTs  10 

GET=  .35  CMS 

QTRlb  < CMS) -  2.55  3.57 

TTRIB  <C)=  9.8  9.8 

ITRIB=  9  17 

KTRIB=  35  10 

QUO  (CMS)s  6.46  11.24 

I *D=  9  10 

K«0  =  28  3S 

CSME=  .26E-05  M/S  SRO=  0.  <C  M*M*M)/<S  M*M>  ET=  -1.9  C 

TO=  -2.2  C  «As  1.6  M/S  PHI=  1.57  RAD 

RATIO  OF  SPACE  TO  TIME  INTEGRATED  VOLUME:  100.000  PER  CENT 

IL=  2  KT=  3  EL=2748.53  M 

22  IM) 

1234  579  12  15 

-1.2285-1.2285-1.2285-1.2284-1.2284-1.2283-1.2283-1.2282-1.2282 


.‘■I 


USER  GUIDE  EXAMPLE  —  DILLON  RESERVOIR.  COLORADO  —  2/18/81 
ELTM=  215.00000  DAYS 


U  <M/S> 


1 

2 

3 

4 

5 

7 

9 

12 

15 

3 

.0007 

-.0037 

-.0102 

-.0200 

-.0255 

-.0173 

-.0179 

-.0192 

-.0082 

4 

.0007 

.0047 

.0014 

-.0052 

-.0098 

-.0031 

-.0012 

.0011 

.0070 

5 

.0007 

.0041 

.0042 

.0023 

-.0004 

.0011 

.0014 

.0054 

.0094 

6 

.0007 

.0011 

.0039 

.0053 

.0043 

.0020 

.0017 

.0061 

.0088 

7 

.  0007 

-.0005 

.0034 

.0062 

.0059 

.0021 

.0017 

.0062 

.0063 

8 

.  0007 

-.0003 

.0034 

.0064 

.0063 

.0021 

.0017 

.0062 

-.0001 

9 

.0007 

.0011 

.0037 

.0065 

.0063 

.0021 

.0017 

.0062 

-.0113 

10 

.0007 

.0033 

.0042 

.0065 

.0063 

.0021 

.0017 

.0062 

-.0217 

11 

.0007 

.0062 

.0049 

.0066 

.0062 

.0021 

.0017 

.0062 

-.0241 

12 

.0007 

.0081 

.0051 

.0066 

.0062 

.0021 

.0017 

.0062 

-.0273 

13 

0.0000 

0.0000 

.0053 

.0  066 

.0062 

.0021 

.0017 

.0063 

-.0296 

14 

0. 0000 

0.0000 

.0057 

.0067 

.0062 

.0021 

.0017 

.0064 

0.0000 

15 

0.0000 

0.0000 

.0062 

.0068 

.0062 

.0021 

.0017 

.0064 

0.0000 

16 

0.0000 

0.0000 

0.0000 

.0072 

.0062 

.0021 

.0017 

.0064 

0.0000 

17 

0.0000 

0.0000 

0.0000 

.0074 

.0062 

.0021 

.0017 

.0063 

0.0000 

18 

0. 0000 

0.0000 

0.0000 

.007A 

.0062 

.0021 

.0017 

.0061 

0.0000 

19 

0.0000 

0.0000 

0.0000 

.0083 

.0063 

.0021 

.0017 

.0060 

0.0000 

20 

0.0000 

0.0000 

0.0000 

0.0000 

.0066 

.0020 

.0017 

.0059 

0.0000 

21 

0.0000 

0.0000 

0.0000 

0.0000 

.0068 

.0020 

.0017 

.0059 

0.0000 

22 

0.0000 

0.0000 

0.0000 

0.0000 

.0087 

.0020 

.0017 

.0058 

0.0000 

23 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

.0020 

.0017 

.0058 

0.0000 

24 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

.0020 

.0017 

.0058 

0.0000 

25 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

.0020 

.0017 

.0057 

0.0000 

26 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

.0020 

.0017 

.0056 

0.0000 

27 

0.0000 

0.0000 

0.0000 

o.o«eo 

0.0000 

.0021 

.0017 

.0056 

0.0000 

28 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

.0020 

.0018 

.0056 

0.0000 

29 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

•  0019 

.0017 

.0055 

0.0000 

30 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

.0019 

.0017 

.0054 

0.0000 

31 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

.0020 

.0017 

.0052 

0.0000 

32 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

.0018 

.0017 

.0044 

0.0000 

33 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

.0021 

0.0000 

0.0000 

34 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

.0061 

0.0000 

0.0000 

35 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

.0226 

0.0000 

0.0000 
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USER  GUIDE  EXAMPLE  —  DILLON  RESERVOIR#  COLORADO  —  2/18/81 
ELTM=  215.00000  CAYS 

U  (MM/S) 


1 

2 

3 

4 

5 

7 

9  12 

15 

3 

0.0000 

.0349 

.0468 

.0617 

.0424 

.0127 

.0143  -.0295 

-.0515 

4 

0.0000 

.0245 

.0617 

.0891 

.0624 

.0023 

.0169  -.0422 

-.0684 

5 

0.0000 

.0121 

.0628 

.0969 

.0720 

“.0040 

.0172  -.0535 

-.0781 

6 

0.0000 

.0126 

.0585 

.0951 

.0756 

'.0041 

.0168  -.0664 

-.0889 

7 

0.0000 

.0280 

.0554 

.0927 

.0778 

-.0019 

.0165  -.0815 

-.1027 

8 

0.0000 

.0530 

.0534 

.0900 

.0793 

.0007 

.0163  -.0946 

-.1141 

9 

0.0000 

.0807 

.0535 

.0891 

.0814 

.0034 

.0161  -.1035 

-.1212 

10 

0.0000 

.1011 

.0575 

.0920 

.0846 

.0060 

.0159  -.1088 

-.1265 

11 

0.0000 

.0889 

.0665 

.0990 

.0887 

.0065 

.0156  -.1129 

-.1167 

12 

0. 0000 

0.0000 

.0843 

.1094 

.0936 

.0108 

.0154  -.1194 

-.0798 

13 

0.0000 

o.cooo 

.0877 

.1247 

.0990 

.0130 

.0152  -.1296 

-.0216 

14 

0.0000 

0.0000 

.0711 

.1432 

.1023 

.0149 

.0149  -.1418 

-.0157 

15 

0.0000 

0.0000 

0.0000 

.1610 

.1016 

.0166 

.0146  -.1529 

-.0089 

16 

0.0000 

0.0000 

0.0000 

.1589 

.1003 

.0182 

.0142  -.1618 

0.0000 

17 

0.0000 

0.0000 

0.0000 

.1489 

.1014 

.0196 

.0139  -.1693 

0.0000 

18 

0. 0000 

0.0000 

0.0000 

.1189 

.1042 

.0208 

.0137  -.1726 

0.0000 

19 

0.0000 

0.0000 

0.0000 

0.0000 

.1096 

.0222 

.0135  -.1704 

0.0000 

20 

0.0000 

0.0000 

0.0000 

0.0000 

.1013 

.0243 

.0135  -.1654 

0.0000 

21 

0.0000 

0.0000 

0.0000 

0.0000 

.0954 

.0274 

.0138  -.1611 

0.0000 

22 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

.0302 

.0146  -.1582 

0.0000 

23 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

.0325 

.0159  -.1518 

0.0000 

24 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

.0330 

.0161  -.1343 

0.0000 

25 

0.0000 

0.0000 

o.cooo 

0.0000 

0.0000 

.0321 

.0163  -.1170 

0.0000 

26 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

.0350 

.0181  -.1072 

0.0000 

27 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

.0420 

.0211  -.0990 

0.0000 

28 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

.0388 

.0083  -.0923 

0.0000 

29 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

.0335 

.0091  -.0872 

0.0000 

30 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

.0259 

.0104  -.0838 

0.0000 

31 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

.0154 

.0130  -.0848 

0.0000 

32 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

.0181  -.0907 

0.0000 

33 

0.0000 

0.0000 

'0.0000 

0.0000 

0.0000 

0.0000 

.0212  -.0842 

0.0000 

34 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

-.0107  -.0637 

0.0000 

C6 


USER  GUIDE  EXAMPLE  — 
ELTMs  215.00000  DAYS 


DILLON  RESERVOIR*  COLORADO 


2/18/81 


T2  <CJ 


1 

2 

3 

4 

5 

7 

9 

12 

15 

3 

0.0 

12.1 

12.0 

12.0 

12.0 

12.0 

11.9 

11.9 

11.7 

4 

0.0 

11.3 

11.3 

11.2 

11.0 

10.9 

10.7 

10.6 

10.6 

5 

0.0 

10.9 

10.8 

10.7 

10.7 

10.6 

10.5 

10.5 

10.5 

6 

0.0 

10.6 

10.6 

10.6 

10.5 

10.5 

10.5 

10.5 

10.5 

7 

0.0 

10.5 

10.5 

10.5 

10.5 

10.5 

10.5 

10.5 

10.5 

6 

0.0 

10.5 

10.5 

10.5 

10.5 

10.5 

10.5 

10.5 

10.5 

9 

0.0 

10.5 

10.5 

10.5 

10.5 

10.5 

10.5 

10.5 

10.5 

10 

0.0 

10.5 

10.5 

10.5 

10.5 

10.5 

10.5 

10.5 

10.5 

11 

0.0 

10.5 

10.5 

10.5 

10.5 

10.5 

10.5 

10.5 

10.5 

12 

0.0 

10.5 

10.5 

10.5 

10.5 

10.5 

10.5 

10.5 

10.5 

13 

0.0 

0.0 

10.5 

10.5 

10.5 

10.5 

10.5 

10.5 

10.5 

14 

0.0 

0.0 

10.5 

10.5 

10.5 

10.5 

10.  5 

10.5 

10.5 

15 

0.0 

0.0 

10.5 

10.5 

10.5 

10.5 

10.5 

10.5 

10.5 

16 

0.0 

0.0 

0.0 

10.5 

10.5 

10.5 

10.5 

10.5 

10.5 

17 

0.0 

0.0 

0.0 

10.5 

10.5 

10.5 

10.5 

10.5 

0.0 

18 

0.0 

0.0 

0.0 

10.5 

10.5 

10.5 

10.5 

10.5 

0.0 

19 

0.0 

0.0 

0.0 

10.5 

10.5 

10.5 

10.5 

10.5 

0.0 

20 

0.0 

0.0 

0.0 

0.0 

10.5 

10.5 

10.5 

10.5 

0.0 

21 

0.0 

0.0 

0.0 

0.0 

10.5 

10.5 

10.5 

10.5 

0.0 

22 

0.0 

0.0 

0.0 

0.0 

10.5 

10.5 

10.5 

10.5 

0.0 

23 

0.0 

0.0 

0.0 

0.0 

0.0 

10.5 

10.5 

10.5 

0.0 

24 

0.0 

0.0 

0.0 

0.0 

0.0 

10.5 

10.5 

10.5 

0.0 

25 

0.0 

0.0 

0.0 

0.0 

0.0 

10.5 

10.5 

10.5 

0.0 

26 

0.0 

0.0 

0.0 

0.0 

0.0 

10.5 

10.5 

10.5 

0.0 

27 

0.0 

0.0 

0.0 

0.0 

0.0 

10.5 

10.5 

10.5 

0.0 

28 

0.0 

0.0 

0.0 

0.0 

0.0 

10.5 

10.5 

10.5 

0.0 

29 

0.0 

0.0 

0.0 

0.0 

0.0 

10.5 

10.5 

10.5 

0.0 

30 

0.0 

0.0 

0.0 

0.0 

0.0 

10.5 

10.5 

10.5 

0.0 

31 

0.0 

0.0 

0.0 

0.0 

0.0 

10.5 

10.5 

10.5 

0.0 

32 

0.0 

0.0 

0.0 

0.0 

0.0 

10.5 

10.5 

10.5 

0.0 

33 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

10.5 

10.5 

0.0 

34 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

10.3 

10.5 

0.0 

35 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

9.9 

10.5 

0.0 
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USER  GUIOE  EXAMPLE  —  DILLON  RESERVOIR,  COLORADO  —  2/18/fll 
ELTMs  215.00000  DAYS 

C2  FOR  CONSTITUENT  1 


1 

2 

3 

4 

5 

7 

9 

12 

15 

3 

0.000 

.020 

.020 

.02  0 

.020 

.020 

.020 

.020 

.020 

4 

0.000 

.021 

.020 

.020 

.020 

.020 

.020 

.020 

.020 

5 

0.000 

.021 

.020 

.020 

.020 

.020 

.020 

.020 

.020 

6 

0.000 

.021 

.020 

.020 

.020 

.020 

.020 

.020 

.020 

7 

0.000 

.021 

.020 

.020 

.020 

.020 

.020 

.020 

.020 

8 

0.000 

.021 

.020 

.020 

.020 

.020 

.020 

.020 

.020 

9 

0.000 

.021 

.020 

.020 

.020 

.020 

.020 

.020 

.020 

10 

0.000 

.021 

.020 

.020 

.020 

.020 

.020 

.020 

.021 

11 

0.000 

.020 

.020 

.020 

.020 

.020 

.020 

.020 

.021 

12 

0.000 

.020 

.p20 

.020 

.020 

.020 

.020 

.020 

.021 

13 

0.000 

0.000 

.020 

.020 

.020 

.020 

.020 

.020 

.021 

14 

0.000 

0.000 

.020 

.020 

.020 

.020 

.020 

.020 

.020 

15 

0.000 

0.000 

.020 

.020 

.020 

.020 

.020 

.020 

.020 

16 

0.000 

0.000 

0.000 

.020 

.020 

.020 

.020 

.020 

.020 

17 

0.000 

0.000 

0.000 

.020 

.020 

.020 

.020 

.020 

0.000 

18 

0.000 

0.000 

0.000 

.020 

.020 

.020 

.020 

.020 

0.000 

19 

0.000 

0.000 

0.000 

.020 

.020 

.020 

.020 

.020 

0.000 

20 

0.000 

0.000 

0.000 

0.000 

.020 

.020 

.020 

.020 

0.000 

21 

0.000 

0.000 

0.000 

0.000 

.020 

.020 

.020 

.020 

0.000 

22 

0.000 

0.000 

0.000 

0.000 

.020 

.020 

.020 

.020 

0.000 

23 

0.000 

0.000 

0.000 

0.000 

0.000 

.020 

.020 

.020 

0.000 

24 

0.000 

0.000 

0.000 

0.000 

0.000 

.020 

.020 

.020 

0.000 

25 

0.000 

0.000 

0.000 

0.000 

0.000 

.020 

.020 

.020 

0.000 

26 

0.000 

0.000 

0.000 

0.000 

0.000 

.020 

.020 

.020 

0.000 

27 

0.000 

0.000 

0.000 

0.000 

0.000 

.020 

.020 

.020 

0.000 

28 

0.000 

0.000 

0.000 

0.000 
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GENERAL 

LARM2  identifies  two  types  of  errors:  (1)  fatal  errors 
which  stop  the  computation;  (2)  non-fatal  errors  which  permit  the 
computation  to  continue.  The  messages  written  on  the  detection 
of  these  errors  are  discussed  below. 

FATAL  INPUT  ERROR 

This  error  occurs  when  the  input  stream  is  missing  a  re¬ 
quired  card  or  the  deck  is  out  of  order  (see  Appendix  A).  The 
code  checks  the  card  identifiers  (field  0)  and  stops  the  compu¬ 
tation  when  an  error  is  detected.  An  input  error  also  occurs 
when  too  few  or  too  many  outlets,  tributaries,  or  withdrawals 
are  specified.  In  both  cases,  LAR1I2  writes  an  error  message, 
repeats  the  offending  card,  and  stops  the  computation. 

NON-FATAL  ERROR:  AZO  EXCEEDS  AZMAX 
and 

NON-FATAL  ERROR:  DZO  EXCEEDS  DZMAX 

The  variables  AZO  and  DZO  are  set  by  the  user  and  represent 
base  values  of  vertical  dispersion  coefficients  for  the  x-  momentum 
and  for  temperature  and  constituent  equations,  respectively.  These 
variables  are  used  with  the  Richardson  number  to  compute  spatially 
and  temporally  varying  dispersion  coefficients.  The  upper  limits 
on  the  latter,  however,  are  values  based  on  the  horizontal  layer 
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thickness,  h,  and  the  time  step  At  (see  Appendix  E).  The  values 
of  AZO  and  DZO  should  always  be  less  than  these  maximum  values. 

This  error  is  not  serious  enough  to  stop  the  computation,  however, 
since  the  code  always  insures  that  the  computed  values  of  AZ  and 
DZ  do  not  exceed  their  computational  limits. 

FATAL  COMPUTATION  ERROR:  VRATIO  EXCEEDS  ONE  PER  CENT  DEVIATION 
LIMIT 

This  message  indicates  one  of  two  conditions.  Most  frequently, 
a  computational  instability  has  occurred  and  the  At  needs  to  be 
decreased. 

On  very  long  simulations  with  short  At's,  this  volume  balance 
variable  may  deviate  somewhat  from  its  perfect  value  of  exactly 
100  due  to  roundoff  error.  This  error  does  not  affect  the  accu¬ 
racy  of  the  results.  In  cases  where  this  roundoff  error  exceeds 
the  one  per  cent  limit  by  gradually  increasing  over  time,  the 
computation  may  be  restarted  part  way  through  the  simulation. 

FATAL  TVDS  ERROR  AT  -  DAYS 

This  error  message  is  coded  in  the  TVDS  subroutine  to  detect 
the  ca^se  where  LARM2  asks  for  time-varying  data  that  is  not  avail¬ 
able  if,  for  example,  a  simulation  were  set  up  for  a  370-day  run, 
with  only  one  year  of  data  available.  The  illegal  time  request 
is  printed  with  the  error  message,  and  the  computation  is  stopped. 


The  glossary  lists  only  those  FORTRAN  variables  found  on 
the  LARM2  printed  output,  in  the  FORTRAN  statements  labelled 
CHANGE,  in  the  TVDS  arguments,  and  in  this  report.  Those 
variables  discussed  in  the  card  description  Appendix  A  are  not 
repeated  here. 

A:  vector  containing  subdiagonal  (backward)  coefficients  for 
the  tridiagonal  algorithm  solutions  of  Z,  T,  and  C 

AR:  waterbody  area  by  layer,  m2  (vector) 

AS:  A  Storage:  array  for  storing  A  vectors  for  use  in  con¬ 
stituent  concentration  computations 

AX:  Ax,  m2/s;  x-direction  momentum  dispersion  coefficient 

AZ:  A  ,  m2/s;  dispersion  coefficient,  z-direction,  x-momentum 

equation  (array) 

AZMAX:  maximum  permitted  value  of  A_  (0,85*i,h2/At) ,  m2/s 

z 

AZMIN:  minimum  value  of  A  (molecular  value),  m2/s  (scalar) 

z 

AZO:  base  value  for  A  computation,  m2/s  (scalar) 

z 

B:  B,  m,  cell  width  (array) 

BETA:  8,  fraction  of  incident  solar  radiation  absorbed  at  the 

water  surface 

BH1:  B’h,  m2 ,  cross-sectional  area  of  a  cel  1  at  the  current  or 
previous  time  step 

BH2:  B-h,  m2 ,  cross-sectional  area  of  a  cell  at  the  previous  or 
current  time  step 

C:  vector  containing  superdiagonal  (forward)  coefficients  for 
the  tridiagonal  algorithm  solution  of  Z  and  T 

CHZY:  Chezy  resistance  coefficient,  mfi/s 
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CIN:  upstream  inflow  constituent  concentration,  mg* l-1  (vector) 

CS:  C  Storage:  array  for  storing  C  vectors  for  use  in  con¬ 

stituent  concentration  computations 

CSHE:  k,  kinematic  coefficient  of  surface  heat  exchange, 
m*s” 

CTRIB:  array  containing  tributary  constituent  concentrations, 

mg* l-1 

CZ:  C*,  wind  resistance  coefficient;  varies  with  wind  speed 

Cl:  array  containing  constituent  concentrations  at  the  current 

or  previous  time  step,  mg-l"1 

C2:  array  containing  constituent  concentrations  at  the  current 

or  previous  time  step,  mg*l_1 

D:  vector  containing  right-hand  side  for  the  tridiagonal 

algorithm  solution  of  Z  and  T 

DX:  D  ,  m2/s;  dispersion  coefficient,  x-direction,  heat  and 
constituent  balance  equation  (array) 

DXI:  base  value  of  Dx,  m2/s  (scalar) 

DZ:  Dz,  m2/s;  dispersion  coefficient,  z-direction,  heat  and 

constituent  balance  equation  (array) 

DZMAX :  maximum  permitted  value  of  D  (0.85-i*h2 /At) ,  m2/s 
(scalar) 

DZMIN:  minimum  value  of  D_  (molecular  value),  m2/s  (scalar) 

z 

DZO:  base  value  of  Dz  computation,  m2/s  (scalar) 

EL:  water  surface  elevation,  m 

ELTM:  elapsed  time,  days;  sum  of  NAt  and  TMSTRT ;  time  since 

beginning  of  calendar  year 

ET:  E,  equilibrium  temperature  of  surface  heat  exchange, °C 
G:  g,  m/s2 ;  gravitational  constant 

GAMMA:  y,  m-1 ;  exponential  decay  constant  for  absorption  of 

solar  radiation  with  depth 
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HIN:  H  Initial,  m;  retains  original  value  of  h  (layer  thickness) 
for  computation  of  varying  h  in  top  layer  of  cells  from  Z 
and  also  used  in  reinitializing  a  layer  of  cells  immedi¬ 
ately  after  a  layer  has  been  added  or  subtracted  due  to 
reservoir  volume  changes 

HN:  H  , °C-m?s_1  or  mg •  l-.1m?s- 1 ;  net  rate  of  heat  or  constituent 

addition;  source  term  for  heat  balance  or  constituent 
equation 

HI:  h,  m;  vertical  slice  thickness  at  the  current  or  previous 

time  step.  (Only  the  top  layer  of  cells  has  varying  thickness 
due  to  elevation  changes,  however,  the  location  of  the  top 
layer  of  cells  may  vary  due  to  extreme  elevation  changes) 

H2 :  h,  m;  vertical  slice  thickness  at  the  current  or  previous 

time  step 

I:  counter  for  sweeping  across  (longitudinally)  the  reservoir 

grid  in  x-direction  DO-loops 

IB:  I  Begin:  variable  locating  the  first  active  cell  in  each 

line  for  use  in  the  tridiagonal  solution  of  the  heat  bal¬ 
ance  equation;  IB  is  found  in  an  initial  program  step  and 
stored  in  the  array  LC 

IE:  I  End:  variable  locating  the  last  active  cell  in  each 

line  for  use  in  the  tridiagonal  solution  of  the  heat  bal¬ 
ance  equation;  IE  is  found  in  an  initial  program  step  and 
stored  in  the  array  LC 

IFLAG:  logic  flag  for  subroutine  TVDS;  if  IFLAG=0 ;  instruct 

TVDS  to  read  time-varying  boundary  condition  data;  if 
IFLAG=1,  instruct  TVDS  to  provide  boundary  condition 
data  for  a  particular  time;  if  IFLAG=-1,  TVDS  error 
and  computation  stopped 

I FORM:  control  variable;  if  IFORM=0,  terminal-oriented  display 

is  produced  (80-character  width);  if  IF0RM=1 ,  full- 
width  (132-character)  display  is  produced 

IL:  I  Left:  I-index  of  current  left-hand  segment 

IMAX:  I  MAXimum:  I-index  of  the  rightmost  segment 

IMAXM1:  I  MAXimum  Minus  1:  I-index  of  the  next-to-rightmost 

segment 

IMAXM2 :  I  MAXimum  Minus  2:  I-index  of  the  second-to-rightmost 

segment 
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ISC:  Integer  Segment  Coordinates:  vector  containing  K-index 

of  bottom  active  cell  of  each  segment 

IY :  Integer  Year:  year  number 

JC:  constituent  counter 

K:  counter  for  sweeping  down  (vertically)  the  reservoir  grid 

in  z-direction  DO-loops 

KB:  K  Bottom:  K-index  of  bottom  active  cell  in  a  particular 

segment;  KB  is  found  in  array  ISC 

KMAX:  K-MAXimum:  K-index  of  the  bottom  layer  of  cells 

KMAXM1:  K-MAXimum  Minus  1:  K-index  of  the  next-to-bottom 

layer  of  cells 

KOUT:  computed  layer  number  at  which  each  tributary  enters 

the  LARM2  grid  (vector) 

KT:  K-Top:  K-index  of  the  currently  active  top  layer  of  cells 

L:  integer  variable  denoting  active  cells  (L-l)  or  inactive 
cells  (L=0) 

LC:  Layer  Coordinates:  array  containing  (1)  layer  number  (K) , 
(2)  IB,  and  (3)  IE  for  each  active  layer  of  cells 

N:  counter  for  time  steps 

P:  pressure,  Pa 

PHI:  4> ,  degrees;  wind  direction 

QET:  Q  Evaporation  Total,  m3/s;  evaporation  rate 

QIN:  upstream  inflow  rate,  m3s_1 

QOUT:  vector  containing  downstream  outflow  rates 

QTRIB:  Q  TRIButary,  m3/s;  vector  containing  tributary  inflow 
rates 

QWD:  Q  WithDrawal,  m3/s;  vector  containing  withdrawal  rates 
RHO:  p,  kg/m3;  density  of  water 

RHOA:  pa,  kg/m3;  density  of  air  (used  in  wind  stress  computa¬ 
tion  and  assumed  constant  at  1.25  kg/m3) 


RR:  constituent  reaction  rate,  s-1 


SRO:  solar  radiation,  (°C-m3/s)/m2 

T:  vector  containing  the  results  of  the  layer-by-layer  solution 

of  the  implicit  temperature  equation 

TAPE5:  read  device 

TAPE6:  write  device 

TAPE61 :  binary  write  device 

TD:  T^^C:  dew  point  temperature 

TIN:  upstream  inflow  temperature , °C 

TMSTRT:  time  of  year  in  Julian  days  when  simulation  starts; 

read  in  as  days,  converted  internally  to  seconds 

TTRIB:  Temperature  of  TRIButary , °C ;  vector  containing  tribu¬ 
tary  temperatures 

T1 :  T,°C;  array  containing  temperatures  at  the  current  or 

previous  time  step 

T2:  T,°C;  array  containing  temperatures  at  the  previous  or 

current  time  step 

U:  U,  m/s;  array  containing  x-direction  velocity  components 

UM:  cell  average  x-direction  velocity,  m  s-1 

V:  vector  containing  main  diagonal  (centered)  coefficients  for 

the  tridiagonal  algorithm  solution  of  Z  and  T 

VOL:  reservoir  volume,  m3 ;  based  on  initial  volume  and  sum  of 

inflows  and  outflows 

VOLE:  cumulative  waterbody  volume  by  layer,  m3  (vector) 

VOLP:  reservoir  volume,  m3 ;  computed  from  geometry  and  water 

surface  elevations 

VRATIO:  ratio  of  VOLP  to  VOL,  % 

VS:  v  Storage:  array  for  storing  V  vectors  for  use  in  con¬ 

stituent  concentration  computations 
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W:  m/s;  array  containing  vertical  velocities 

WA:  w  .  m/s;  wind  speed 

a 

WKT:  a  vector  containing  a  corrected  vertical  velocity  for 

the  top  layer  of  cells  that  permits  a  perfect  water 
balance  in  that  layer 

WM:  cell  average  z-direction  velocity  m  s’*1 

Zl:  Z,  m;  vector  containing  deviation  of  water  surface  from 

datum  (top  of  active  layer  of  cells)  for  each  column; 
positive  downwards.  The"l"  denotes  a  particular  time 
level 

Z2:  Z,  m;  vector  containing  deviation  of  water  surface  from 
datum  (top  of  active  layer  of  cells)  for  each  column; 
positive  downwards.  The "2"  denotes  a  particular  time 
level 
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GIN  is  a  preprocessor  code  that  prepares  bathymetry  (BA) 
cards  for  use  by  LARM2.  GIN  uses  GEDA  output  as  input  and  may 
be  used  to  preview  the  grid  geometry  prior  to  creating  the  BA 
cards.  The  code  is  short  and  is  documented  internally  with 
comment  cards.  Following  are  the  most  important  features  of 
GIN. 


‘  The  basis  grid  parameters,  Ax  and  h,  are  selected 
during  the  GEDA  run.  Output  from  the  final  GEDA  run 
is  modified  by  eliminating  all  output  lines  preceding 
"GEOMETRIC  MODEL  FOR  UNSTEADY  FLOW  PROGRAMS."  This 
modified  GEDA  output  is  the  only  input  to  GIN. 

’  The  variable  IFORM  can  be  set  to  0  or  1  in  line  27 
of  GIN.  If  IFORM=0,  the  grid  geometry  is  previewed 
in  a  form  similar  to  that  presented  by  LARM2.  If 
IF0RM=1,  the  BA  cards  are  written.  The  user  will 
generally  run  GIN  with  IFORM=0  until  a  satisfactory 
grid  is  established,  then  with  IF0RM=1  as  a  final 
step  before  working  with  LARM2. 

*  The  segment  order  may  be  reversed  by  changing  the 
variable  IMAGE. 

*  GIN  may  be  used  to  modify  a  grid  geometry  by 
changing  some  lines  of  the  code.  The  example  for 
Dillon  Reservoir  shows  some  lines  marked  "CHANGE" 
in  columns  73-80  in  DO-loop  170  that  modify  widths 
to  fit  better  the  given  elevation-area-volume  table. 
The  best  compromise  in  fitting  such  tables  may  be 
to  match  areas  in  the  surface  layers  and  volumes 
elsewhere. 

*  The  user  is  responsible  for  obtaining  a  legal  geom¬ 
etry  as  outlined  in  Section  4.2  of  this  report. 

'  GIN  uses  GEDA  output  in  feet  and  converts  all  dimen¬ 
sions  to  metres. 


W  is  a  postprocessor  code  that  uses  the  velocity  fields 
generated  by  LARM2  to  produce  vector  plots  like  that  shown  on 
page  G3.  W  uses  DISSPLA,  a  proprietary  software  package 
created  by  Integrated  Software  Systems  Corporation  of  San  Diego 
and  available  on  many  time-sharing  systems.  W  uses  a  file 
created  by  LARM2  when  the  IPLOT  parameter  on  the  PL  card  is 
set  to  one.  The  code  is  short  and  documented  internally  with 
comment  cards.  However,  familiarity  with  DISSPLA  subroutines 
is  required  to  use  W.  Following  are  some  of  the  most  important 
features  of  W. 


'  LARM2  writes  the  contents  of  the  TITLE  and  GE  cards  to 
TAPE61  in  binary  prior  to  writing  the  coordinates  of 
all  the  LARM2  cells  and  the  velocity  fields.  W  reads 
this  information  on  TAPE51. 

'  W  computes  all  the  scaling  parameters  and  axis  labels 
from  information  written  by  LARM2,  The  user  can  change 
the  axis  sizes  and  location  of  the  plot  on  the  paper  by 
changing  the  values  of  XAXIS,  ZAXIS,  XSTART,  and 
ZSTART  (all  in  inches). 

'  The  velocity  components  are  translated  to  a  displace¬ 
ment  from  the  cell  center  by  multiplying  each  compo¬ 
nent  by  DLT,  a  time  in  seconds.  Different  applications 
yield  different  velocity  magnitudes,  so  the  user  may 
need  to  change  the  value  of  DLT.  It  is  converted  to 
days  and  identified  as  the  scaling  parameter  on  each 
plot . 

'  The  velocity  components  used  by  W  are  averaged  to  the 
cell  center  so  that  U(I,K)  =  (U(I-1,K)  +  U(I,K))/2  and 
f  =  (W( I , K-l)  +  W( I , K) )/2 .  This  averaging  process  may 
result  in  stepped  velocity  profiles  adjacent  to  vt  ti- 
cal  boundaries. 
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The  user  may  select  the  number  of  velocity  plots  by 
manipulating  the  PSTART  and  PEND  parameters. 

W  also  writes  all  the  plot  definition  variables  and 
title  information  to  TAPE6  as  a  debugging  aid. 
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Hydrodynamics  and  Transport  of  Chlorine  in  Panther  Branch  Arm, 
Squaw  Creek  Reservoir  for  Commanche  Peake  S.E.S.,  John  Eric 
Edinger  and  Edward  M.  Buchak,  prepared  for  Texas  Utilities 
Services,  Inc.,  Dallas,  Tex.,  June  1978. 

A  Hydrodynamic.  Two-Dimensional  Reservoir  Model:  Development  and 
Test  Application  to  Sutton  Reservoir  Elk  River.  West  Virginia 
John  Eric  Edinger  and  Edward  M.  Buchak,  prepared  for  U.  S.  Army 
Engineer  Division,  Ohio  River,  Cincinnati,  June  1979. 

Temperature  and  Hydrodynamic  Predictions  for  Center  Hill  Lake 
Using  the  LARM  Two-Dimensional  Computer  Model ,  John  A.  Gordon, 
prepared  for  U.  S.  Army  Engineer  District,  Nashville,  June  1979. 

A  Review  of  Numerical  Reservoir  Hydrodynamic  Modeling,  Billy  H. 
Johnson,  Technical  Report  E-81-2,  U.  S.  Army  Engineer  Waterways 
Experiment  Station,  Vicksburg,  Miss.,  November  1981. 

Hydrothermal  Simulations  of  Commanche  Peake  Safe  Shutdown 

Impoundment ,  Edward  M.  Buchak  and  John  Eric  Edinger,  prepared 
for  Texas  Utilities  Services,  Inc.,  Dallas,  Tex.,  May  1980. 

"An  Evaluation  of  the  LARM  Two-Dimensional  Model  for  Water 
Quality  Management  Purposes,"  John  A.  Gordon,  published  in 
Proceedings  of  the  Symposium  on  Surface-Water  Impoundments 
held  in  June  I960,  Minneapolis,  Minn. ,  by  the  American  Society 
of  Civil  Engineers. 

Thermal  Demonstration  Pursuant  to  Illinois  Pollution  Control 
Board  Rules  and  Regulations  Chapter  3,  Rule  203(i)(10),  Energy 
Impact  Associates,  prepared  for  Illinois  Power  Company, 

Decator,  Ill.,  July  1980. 

Development  Document  Relating  to  the  1979  Temperature  Simulation 
Studies  Battle  River  Reservoir  Phaso  1,  MacLaren  Engineers 
Planners  and  Scientists,  Inc.,  and  J.E.  Edinger  Associates, 

Inc.,  prepared  for  Alberta  Power  Limited,  Edmonton,  Alberta, 
Canada,  October  1980. 

Estuarine  Laterally  Averaged  Numerical  Dynamics:  The  Development 
and  Testing  of  Estuarine  Boundary  Conditions  in  the  LARM  Code, 
John  Eric  Edinger  and  Edward  M.  Buchak,  prepared  for  U.  S.  Army 
Engineer  District,  Savannah,  July  1980. 

Un-ionizad  Ammonia  Distribution  in  Blue  River  Arm  of  Dillon 
Reservoir,  J.E.  Edinger  Associates,  Inc. ,  and  Black  &  Veatch , 
Kansas  City,  Kans. ,  March  1981. 
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In  accordance  with  letter  from  DAEN-RDC,  DAEN-ASI  dated 
22  July  1977,  Subject:  Facsimile  Catalog  Cards  for 
Laboratory  Technical  Publications,  a  facsimile  catalog 
card  in  Library  of  Congress  MARC  format  is  reproduced 
below. 


Buchak,  Edward  M. 

User  guide  for  LARM2  :  A  longitudinal -vertical, 
time-varying  hydrodynamic  reservoir  model  /  by  Edward 
M.  Buchak  and  John  E.  Edinger  (J.E.  Edinger  Associates, 
Inc.)  —  Vicksburg,  Miss.  :  U.S.  Army  Engineer 
Waterways  Experiment  Station  ;  Springfield,  Va.  ; 
available  from  NTIS,  1982. 

91  p.  in  various  pagings  :  ill.  ;  27  cm.  — 
(Instruction  report  ;  E-82-3) 

Cover  title. 

"December  1982." 

Final  report. 

"Prepared  for  Office,  Chief  of  Engineers,  U.S.  Army 
under  Contract  No.  DACW39-78-C-0047  (Work  Jnit  31593; 
Task  IA.4) ." 

"Monitored  by  Hydraulics  Laboratory,  U.S.  Army 
Engineer  Waterways  Experiment  Station." 

At  head  of  title:  Environmental  6  Water  Quality 
Operational  Studies. 

Bibliography:  p.  32. 


Buchak,  Edward  M. 

User  guide  for  LARM2  :  A  longitudinal-vertical  :  ...  1982 

(Card  2) 

1.  Computer  programs.  2.  Hydrodynamics. 

3.  LARM2  (Computer  program).  4.  Mathematical  models. 

5.  Reservoirs.  6.  Water  quality.  I.  Edinger, 

John  E.  II.  United  States.  Army.  Corps  of  Engineers. 

Office  of  the  Chief  of  Engineers.  III.  Environmental 
8  Water  Quality  Operational  Studies.  IV.  J.E.  Edinger 
Associates,  Inc.  V.  U.S.  Army  Engineer  Waterways 
Experiment  Station.  Hydraulics  Laboratory.  VI.  Title 
VII.  Series:  Instruction  report  (U.S.  Army  Engineer 
Waterways  Experiment  Station)  ;  E-82-3. 

TA7.W34i  no. E-82-3 


